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niversality in the effective pair interaction of d-shell 
etals - Compressibility and vacancy formation energy 
' 3J-liquid metals 
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stract. Like liquid alkali metals d-shcll liquid metals show scaling behaviour in structure and 
iraction potentials. A realistic interaction potential model, properly parametrized can reasonably 
icribe the universality in the isothermal compressibility and vacancy formation energy of 3d- 
aid metals in electron ion plasma model. 

ywords. Universality; effective pair-potential; 3rf-liquid metals. 
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Introduction 

is now well established that for liquid alkali metals, the structure factors and pair 
itribution functions when suitably expressed scaling behaviour follows very well [1]. In 
nsition metals, in spite of the fact that the incomplete d-shell contributes substantially 
the properties, the liquid structure factors can be put into the same class of simple 
uid metals such as Al [2]. A large number of diffraction data on the liquid structure 
:tors of these metals are now available in the literature [2]. A plot of structure factors, 
2) versus <2/<2max where 'max' refers to the value of wave vector at which S(Q] is 
iximum is shown in figure 1 and it is evident that the structure factors scale reasonably 
ill. It is therefore expected that the basic function i.e., the intermolecular potential 
srgy function should scale similarly and thereby show certain universality in the 
srage behaviour of different properties. But unfortunately, all the experimental 
ucture factor data, accessible from diffraction experiments, are devoid of accurate low 
data - in fact the data below Q ~ 1.0 A" 1 is rarely available. Since the accuracy of the 
'ective pair potential extracted from diffraction data depends strongly on the low Q 
ta [3], it is not easy to use the integral equations to extract the effective pair potential 
d see how the scaling property in the potential operates in these liquids. We have 
irefore addressed the problem somewhat indirectly. Earlier works from this laboratory 
ve shown that at least for the liquid 3d-transition metals, the simplified Wills-Harrison 
)del potential [4] in slightly modified form [5, 6] works quite well as far as the structure 
d thermodynamic properties are concerned. We have therefore thought it useful to test 
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Figure 1. Scaling of experimental structure factors [2] of some SaMiquid transition 
metal (all at or near their melting points). 



members of 3fi?-series too scale properly. We have then obtained the structure factors 
using these potentials following the method described in ref. [6] and then Fourier trans- 
formed to obtain pair distribution functions. Both structure factors and pair distribution 
functions agree well with experimental data. It is also tested how these computed 
structure factors and pair distribution functions scale. In 3, the scaling behaviours via 
the isothermal compressibility, a quantity directly linked to the long wave limit of 
structure factor are also investigated. 

In 4 the vacancy formation energies of several members of the 3d-metal series are 
studied via Minchin et al [7] formalism. These authors have given the Faber's formula [8] 
for the vacancy formation energies in a direct space formalism, and have applied to 
mono-, di- and polyvalent metals. More recently Bhatia and March [9] have used the 
Minchin et al formalism to estimate this quantity for solid argon at the melting tempera- 
ture. Das et al [10] successfully applied the method for liquid alkali and several other 
metals at or near their m.p.'s. We have reinvestigated this formalism for several liquid 
metals of 3d-transitibn series using both one component plasma (OCP) and charged hard 
sphere (CHS) plasma models for reference. It is evident that vacancy formation energies 
of several members of 3d-series too scale though in somewhat limited sense. 



2. W-H model for effective potential of M-liauid metals and liouid structure 
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Liquid metals 
with solid 
state structure 


N c 


z d 


ra 
(au) 


r c 
(au) 


Ti (hep) 
Mn (bcc) 
Fe(bcc) 
Co (hep) 
Ni (fee) 
Cu (fee) 


12 
8 
8 
12 
12 
12 


2.5 
5.5 
6.5 
7.5 
8.5 
9.5 


2.04 
1.63 
1.51 
1.44 
1.34 
1.266 


1.88 
2.00 
1.61 
1.48 
1.36 
1.09 



(b) 



Liquid 
metals 


Packing 
density 
r) 


Charge 
parameter 

r 


Ti 
Mn 
Fe 
Co 
Ni 
Cu 


0.437 
0.449 
0.438 
0.444 
0.442 
0.450 


110 
110 
110 
110 
110 
120 



repulsive contribution u dd (r], from the shift of the d-band centres 
(r) = u s (r) + u d (r] + u dd (r}. 



(1) 



The s-electron contribution to the pair potential is given by the Ashcroft empty core form 
of the bare pseudo-potential and in terms of number of s-electrons, and Ashcroft core 
radius, r c the detailed expression is available in ref. [6]. In the construction of dielectric 
screening function same Geldart-Vosko (GV) [11] local field factor as in ref. [6] was 
used. The ^/-electron contributions are expressed in terms of the number of ^-electrons, 
Zj, and the d-state radii r d \ N c is the nearest neighbour coordination number and 
expressions are also given in refs [4] and [6]. 

For all the 3d-transition metals Z s = 1 .5 is used in the calculation. Zd and r d are 
obtained as in ref. [4]. Values of r c derived from the free atom s-state ionization potential 
lie between the two sets of values given by Wills and Harrison as pointed out earlier [6] 
and same values are used in the present calculations. The input parameters are listed in 
table 1. 

The computed effective pair potentials of several liquid 3^-metals are shown in 
figure 2(a). The depth and position of the potential minimum are affected by the choice 
of r c s, and near the minimum the potentials are dominated by the ^-electron contribution. 
The potentials in general are stronger than typical liquid alkali metal potentials but the 
softness of the core potential appears almost alike. The Friedel oscillations [6] are almost 
absent. For comparison we have also shown the effective pair potential of liquid sodium 
at the melting point (m.p.) [12]. In figure 2(b) the scaled property of the effective pair 
potential (in respect of position and depth) for a few members of 3d-metal series is 
shown. The potentials are reasonably scaled. 

The structure factor S(<2)'s obtained from effective pair potential through low order 
perturbation theory like random phase approximation (RPA) agree extremely well with 
experimental data [6]. When Fourier transformed we get pair distribution function, g(r) 
which too agrees well with experimental data. Both S(Q] and g(r) are seen to scale 
surprisingly well. The scaling behaviour of S(Q] and g(r) in respect of the first peak 
height are depicted in figures 3(a) and (b) respectively. These functions scale reasonably 

at hicrhpr tp.mnf.ratiirps fficrnrp. 
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Figure 2. (a) Effective pair-potentials (in au) of six 3d-liquid rnetals alongwith that 
of Na (r c = 1 .69 au and with GV local field); (b) Scaled effective pair-potential of six 
3c?-liquid metals. 
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Figure 3. (a) Calculated S(Q] vs. Q/Q m of six 3d-Hquid metals at T/T mp = 1.0, 
1.11, 1.25. At higher temperatures S(Q] is given for three metals e.g., Fe, Co, Ni. (b) 
Calculated g(r) vs. r/r m of six 3rf-liquid metals (at or near their melting points). 



3. Isothermal compressibility 

Isothermal compressibility, XT is related with the structure factor as [3] 
5(0) = 



(2) 

where 5(0) is the long wave limit of the structure factor 5(fi). Using the compressibility 
equation of state the scaling behaviour of the compressibility XT, of the 3d-transition 
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0.443 
0.433 
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0.467 


0.463 
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0.468 
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liquid metals can be described in the form 

/oo 

5(0) = 1+ 4?rp / [g(r, T] -<l]r 2 dr = 1 + h(Q) (3a) 

Jo 



= 24V2 [g(r*,r)-l]r* 2 dr*, (3b) 



or 

h(0} 



which is a constant at or near m.p. because it depends only on the pair distribution 
function g(r*, T} which is scaled well. In eq. 3(b) r* = r/r m , r m being the position of first 
g(r} maximum and 77 is the hard core packing density. Further eq. 3(a) can be put into the 
form 

/oo 

5(0} = 1 + 24A/277 mp (p/p mp ) / [g(r*, T) - l]r* 2 dr*. (3c) 

Jo 

This equation indicates that variation of 5(0} as a function of p/p mp at higher 
temperatures should be similar for different members of the series. 

Now the long wave limit structure factor 5(0} using the RPA expression [13] for 
transition metals requires the inclusion of ^-electron contribution [6]. Thus, 

5-^0)= S ~ HO) +^0)4-^(0), (4) 

where ^(O) is the contribution for the reference system and in the present case it is 
charged hard sphere (CHS) plasma which depends on parameters 77 and F [14]. 77 and F 
are taken from ref. [6]. S~ l (0}, S^ ] (0} are the contributions for s- and ^-electrons 
respectively in the long wave limit [6]. 

The calculated scaled values of 5(0} are depicted in table 2. The experimental values 
based on Waseda's work [2] are also shown. We also give the values for the liquid alkali 
metals in charged hard sphere plasma model [1]. It is to be noted that all the transition 
metals are scaled very well in long wave limit. The average value for 3d-transition series 
are slightly smaller in magnitude in comparison with the liquid alkali metals. The 
experimental results also show similar difference between liquid transition and alkali 
metals. In figure 4 we have shown variation of S(0)/S(0) mp as a function of p/p m p, where 
n refers to number densitv at elevated temoerature and m.o.. the meltins ooint. F's at 
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Figure 4. Plot of S(0)/S(0) mp vs. p/p mp ; Ti (o o), Mn (D D), Fe (A A), Co (+ +), 
Ni (x x), Cu (). 



higher temperatures are obtained by assuming constant Z s and r/s from an expression 
given by Waseda [2]. For the temperature range considered (7"/r mp = 1.25 maximum), 
the variation of 5(0) is linear and almost similar for different members of 3d-series except 
liquid copper. On various aspects copper appears to behave somewhat differently in 
comparison with other members [4]. The magnitude of the potential minimum is lowest 
for copper and it is lower than average thermal energy, k%T at m.p. 

4. Vacancy formation energy 

A method for the vacancy formation energy of liquid metals involving pseudopotential 
and total pressure was evolved from this laboratory [10]. The expression could be 
conveniently used to compute vacancy formation energy of several simple liquid metals. 
The method here is extended to liquid 3d-metals. Using the direct space formulation of 
Minchin et al [7] (which neglects atomic relaxation around a vacant lattice site) and total 
pressure of the liquid metals, the expression of vacancy formation energy, E v is 
conveniently written as 



-p(du( P }/dp}-(p/2}J M (r,, 
- (p 2 /2] f(d(u(r,p}g(r)}/dp}dr+(p 2 /2} f (dg(r)/dp}u(r,p}dr r (5) 

/ / 

where the terms have their usual meaning as in ref. [10]. Following Bhatia and March [9] 
we use the Ornstein-Zernike (OZ) relation for the direct correlation function (DCF) and 
after some manipulation we get the basic equation 

; = -!- ( P /k E T}(du(p}/dp} - (p/2)[(dS- l (0)/dp) + (dc(r = 0}/dp}} 
( P 2 /2k B T] f(dg(r)/dp}u(r,p)dr, (6) 

J 

where c(r = 0) is the direct correlation function evaluated at r = 0. c(r = 0) can easily 
be evaluated knowing co(r), the reference part and CRp A (r = 0) from the knowledge of 
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Ni 
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9.30 


6.52 
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14.42 


9.53 


8.00 


11.26 


Average values: 






10.09 


7.02 


9.24 


9.64 



indirect s- and ^-electron mediated part in RPA. We use co(r) for two closely related 
electron-ion plasma models for liquid metals, namely, the one component plasma (OCP) 
model of Baus and Hansen [15] as modified by Evans and Sluckin [13] and the charged 
hard sphere (CHS) plasma model [14]. 

The expressions for CQ(T = 0) in the two electron-ion plasma models are given in [10]. 
The CRPA^ = 0) term is given by 

= 0) = -/3/(27r 2 ) j H<2)Q 2 dQ, (7) 

where v(Q) has now been arising from the contributions of s- and ^-electron interactions 
( 2). The expressions for 5o(0) for the two models given in ref. [10] plus the ^-electron 
terms of 3 constitute 5(0). The volume dependent but structure independent energy has 
the same functional form as in ref. [10] but the expression for electron gas free energy, 
F eg is now given as 



z, 
: Y 



0.982/rJ-0.712/r,+0.0311nr s -0.1 15- 



Of) 



y 



(8) 



where the* terms have their usual meaning as in ref. [6]. 

Now in (6) the last term is structure dependent and can be called structural term while 
the rest of the terms are in fact thermodynamic terms [16]. So we can put 



E v /k E T = E*/k E T + EfVfcT 1 . (9) 



It is in general found that the structural term has minor contribution 
B y < i) to the vacancy formation energy and the major part is contributed by 
the thermodynamical terms [16]. Assuming Egelstaff 's uniform fluid model [17] as in ref. 
[10] this can be shown to be approximately the case. We have therefore computed here 
the thermodynamic part only and make the comparison with experimental results. The 
computed results alongwith the experimental values are shown in table 3. The input 
parameters are all taken as in 2 and 3. The results for CHS plasma model are in general 
a little bigger than the values for OCP model. The model results agree fairly well with 

Pramana - J. Phys., Vol. 50, No. 1, January 1998 7 



come out to be about 10. 1 and 7.0 respectively. The results are better than those of Alonso 
and March [18] for liquid transition metals at their m.p's. The average experimental value 
(for 5 metals for which vacancy formation energies are available) is about 9.64. Our 
earlier calculation showed for alkalis (bcc) except for Li an average, E v /k R T as 15.6 and 
12.2 for CHS and OCP models respectively. The average experimental value is 12.2. For 
close packed lattice (both non-transition and transition metals) observed E v /k<B,T ~ 8-9. 
Rashid and March [16] have sought to connect this term i.e., E v /k#T at m.p. to the 
corresponding crystal lattice structure. Thus, the present calculation shows some 
agreement with the proposition of Rashid and March. The graph E v vs T is approximately 
linear (not shown) indicating a moderate scaling behaviour for vacancy formation energy 
for 3d-liquid metals. 



5. Conclusion 

The W-H model for liquid transition metals is easily applicable effective pair potential 
model for the study of various structural and thermodynamic properties of liquid 3d- 
transition metals. The effective pair potentials of this series are scaled well. The structure 
factors 5(2) calculated in RPA and pair distribution functions g(r) also scale surprisingly 
well. The long wave limit of structure factor 5(0) and its equivalent expression r)h(Q}~ is 
found to be constant which indicates again a good scaling behaviour and difference with 
liquid alkali metals is evident. 

The experimental liquid structure factor data for several members of 3<af-series 
depicted in figure 1 however show limited scaling. It is clear that there exist some 
limitations in spite of the fact that all the data shown here are due to Waseda's single-set 
measurements and accuracies in 5(2) and g(r) were claimed about 4%. The coordination 
numbers evaluated from these data varied between 10.6 and 11.6 for different members 
of the series indicating good similarity of nearest neighbour structures in the liquid 
state. The solid state structures, however are different either fcc/hcp or bcc for different 
members of Id-s&ries. Similar is the case with evaluated hard core packing densities 
which varied between 0.44 and 0.45 except for Ti and V indicating once again the 
similarity with limited variations. As pointed out by Waseda and Tamaki [2] the 
measurements in this kind of work, particularly in the data collection, were affected by 
the experimental difficulties arising from the high temperatures due to active chemical 
reactivity of the liquid samples compared to measurements at lower temperature (below 
1373K). In this connection we would also like to mention that the experimental liquid 
structure factor data of alkali metals too show the scaling property with limitations among 
various members of the series [1]. So the scaling depicted in figure 1 for liquid 3d-series 
is not unreasonable. 

We have also applied W-H model potential to obtain vacancy formation energies 
of several liquid 3d-metals in Minchin et al formalism. We have evaluated the 
thermodynamic part only and the average value agrees fairly well with the known 
average values for the 3i/-series. Here also scaling is evident and difference from 
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carriers are found to exist in the (Ba/Sr)-O planes and are coupled to B lg vibrations of the apical 
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1. Introduction 

The nature of the charge carriers in high T c superconductors has been a subject of 
extensive investigations [1^4]. We have proposed that the normal state properties of 
copper-oxide superconductors can be understood on the polaronic model of the charge 
carrier [5]. Recently Zhao et al [2] through the study of isotope effect of the number 
density and London penetration depth have found evidence of polaronic charge carriers in 
L^-jcSr-jcCuCU. Shen and Schrieffer [3] have observed two types of excitations in their 
angle resolved photoemission spectroscopy (ARPES) of Bi(2212), one has a nearly 
Lorentzian shape and is attributed to a quasi-particle and the other has a broad feature 
with an edge type opening near the Fermi surface which they attribute to a collective 
mode of excitation. These line shapes of the excitations are dependent on the temperature, 
level of doping, and momentum of the excitations. Campuzano et al [6] have carried out 
ARPES studies on YBa2Cu30e.9 and have mapped the Fermi surfaces which are in 
agreement with the calculations made using the local-density approximation. Similar 
measurements have been made for Bi (2212) by Marshall et al [7] and they find hole 
pockets centred at (TT, TT). The underdoped sample with less hole has Fermi level crossings 
over a small portion compared to the overdoped sample and its Fermi surface does not 
extend to the Brillioun zone (BZ) edge. The lowest binding energy edge is shifted from 
the Fermi surface by about 20-30 meV as in ref. [3]. There exist several femtosecond 



the nonsuperconductmg YBCO sample. 

We have shown that phenomenological description of the normal state of Cu-O high 
temperature superconductors in terms of q independent polarizability proposed by Varma 
et al [1] follows from the short range cooperative electron-phonon coupling which 
produces correlated polarons in mixed valence compounds [5]. We use this approach and 
show that the ARPES and optical absorption results are satisfactorily explained on the 
basis that both large and small polarons coexist in high-T^ materials and their weight is 
dependent on the doping concentration, temperature and momentum of the excitation. 

2. The model 

We begin with the Holstein molecular crystal model [10] and consider diagonal electron- 
phonon coupling. Taking independent Einstein oscillators having frequency UQ, 



i-a t-ja i 

- A E c ,> C I> M I- + M / 2 E ("? + w o w ?) > (!) 

i-a i 

where fy denotes the bare electron hopping integral from site z to all j sites that are 
available for hopping, A is the electron-lattice coupling constant, HI = cfci, M is the mass 
of the ions and u, their displacement from equilibrium, c(c,- CT ) denotes the creation 
(annhilation) operator for the electron with spin a at site z. 
Using the canonical polaron transformation 

H > H' = &' He ' , (2) 

^A^-^V,,,.^ (3) 



where a* (,-) denotes the creation annihilation operator for the local Einstein oscillators and 

fl i h + fl ' f A \ 

HI = j== . (4) 

^2Mu)Q/h 
We obtain 

H = 



where 

//- = e'V, (6) 

(7) 



is the polaron binding energy and 2e p is the polaron-polaron attraction energy. 
readily shown that there are two different types of solutions, one band type and 
localized. 
The energy of the band type polaron is 



(8) 

It can be 
the other 



and that of the localized state is 

EI = -2J - (A 2 /2Mw 2 ) = -2J - e p . 

We consider that the localized electrons and holes are symmetric about the 
potential, //, so the energy diagram as shown in figure 1 holds. In eqs (9) and 



where (n q ) denotes the set of oscillator occupation numbers, 



r = 



n = 



- I]' 1 , 



(9) 

(10) 

chemical 
(10) 

(11) 

(12) 

(13) 
(14) 



From figure 1 the localized hole states are 2e p below the bottom of the polaron band while 
localized electron states are at the bottom of the band. The hole and electron states arise 
in homogenous mixed valence compound with the metallic ions in A" + and A n+l and 
anions in B m ~ and fi(' M+1 )- states. In a compound like Bi(2212) with some O" 1 ions 
present in the lattice, 2D primitive cell of SrO, assuming that mixed valence cation and 
anion states exist in this plane, will be constituted by two blocks as shown in figure 2. 
Coupling with LO phonon modes with ion displacements normal to the plane along the 
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Figure 2. Location of charges on Sr-O4 plane. 

Bi 3+ -O4 r Cu 2+ bond, will couple 1~, 2~ and 2 + , 1 + states through electron and hole 
transport as discussed in ref. [11] and will lower the energy of the system due to 
asymmetry of the charges in the bond. In this process g is nearly normal to the plane 
containing the X/ and X i+j vectors and leads to near vanishing of r q in eq. (13) so there are 
no Frank-Condon transitions. Consequently there is no narrowing of the polaron 
bandwidth with temperature and no increase in the effective mass of the electron [11]. 
This case was discussed by us in terms of two degenerate states ij)\ = ^/N ex.p(i(f)i) and 
^ I]; = T/N e.xp(i<pu) related by the phase difference (j>\ </>n = 7r/2 and coupled via the 
phonon energy HUJQ. For this case there is no activation energy in the expression for dc 
conductivity. 



3. Band and hopping conductivity 

We now obtain the expressions for optical absorption coefficeint r(w) on the polaronic 
model to explain the results of the ARPES and femtosecond optical absorption studies. 
We assume that the small polaron band states behave approximately the same way as free 
states of large polaron in the continuum approximation. The treatment of Kartheuser [12] 
which considers the electrons initially at rest is applicable to the present case and the 
absorption constant is given by 



1 



2 
3 



} 

m) 






N] 



m Q 
m 12 



(15) 



Here eo is the dielectric constant of vacuum, n(u)} is the index of refraction, c is the 
velocity of light, UJQ is the LO phonon frequency to which the electron is coupled and a is 
the dimensionless coupling constant 

1 e 2 1 /I 1\ 

a = 2~^ OTi^l: -) ( 16 ) 



(17) 



A general approach using density matrix gives 

_. , 1 



sue i -\-j is given oy 

2f,, r\-2(n r^faPM e*p(-g7frft) / IR N 

cr(o;,r) = g*(Q,r) -- -j^, (18) 

(JM3/2) [i + ( 2 ^ Tf) 2 ] 
where 



(19) 

(20) 

Pi = 2(^o6-T 2 ), (20a) 



(i +pf) 2 } 1/2 - {[i +rf] 1/2 - 1}, (20) 

Pi Pi 



and 



f ' l j 



The contribution to dc conductivity is given by 

- ,-1/2 _2_ 



(Ki - ^ +J )(x, - x t+J ) 



"\^>*J 4 ^2 \i i+jJV-i ^i+jJ 

x sech((e,- + ; - ju)/3/2 

x /3-Tj exp(-2G;(#) sin 2 i(0(jc,- - *,-+_/)) tanh ft wo/3/2) . (22) 

Here <:/,, c e are the probability of occupation of sites i and i+jf by hole and electron 
respectively and e,-+ y - /J.= +e p and e/+y /A = e^ as shown in figure 1. In the present 
case with q of the LO phonon normal to the vector (Xi Xf+j), r^ 1 = OJQ and / ~ftwo, we 
obtain for a square lattice with coordination number z = 4, the expression for the dc 
conductivity for the correlated polarons, 

ffc (O, T] = ^-nc e c h a 2 e 2 ujQ/3 sech 2 (e p /3/2), (23) 

where n is the number of ions in the mixed valence state. The optical conductivity can be 
expressed as [1] 

/ rri\ / *TT?\ i / np\ f^)A\ 

The low frequency optical conductivity is obtained from eq. (18) when the electron 
absorbs the photon and emits it without Frank-Condon (FC) transition, 



and e i+ j = e p . (25) 

So f = 0. Then from eq. (19) 

COSh( p /3/2) 



Pramana - J. Phys., Vol. 50, No. 1, January 1998 15 



cr(o,T) 
<r(0,T) 



,0.6 



0-4 



0.0 



p= 8meV 



40 80 120 160 200 
w (meV ) - 



Figure 3. Drude type contribution to optical conductivity from eq. (26) for 
e n = 8meVandr=100K. 



This is Drude type contribution and as shown in figure 3 monotonically decreases with 
increase in u>. 
The second component is obtained for those processes in which the FC transition 

occurs, 

=huj-e e i = e, & = 2e, (27) 



exp(-4e 2 r /2 7?/ft 2 ) 



for u > 2e, 



= for uj < 2e. 



p . 



For u > 2e p , we take 



This gives 
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(28) 
(29) 
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(30) 



The relative contributions to optical conductivity from cr cl (u, T) and o- c2 (a;, T) have to 
be determined from experiment. The nature of cr c2 (a;,r) curve in eq. (30) agrees well 
with that found emperically by Varma et al [1] for 



We now examine the data on optical absorption, photoemission spectroscopy and 
femtosecond absorption in high T c superconductors and show that these can be satisfac- 
torily explained on the model developed above. 

4.1 Optical conductivity 

The optical conductivity of all oxide superconductors seems to be composed of Drude- 
like free carrier peak centered at zero frequency and a broad mid-infrared absorption band 
centered nearly at 2000cm" 1 . Reedyk et al [14] have obtained the optical conductivity of 
Bi(2212) from the reflectance data. The normal state conductivity is dominated by a peak 
at zero frequency representing the Drude contribution. Also there are two broad peaks, 
one near 412cm" 1 and another at 996cm" 1 . As discussed by Cardona et al [15], the 
polarized Raman spectra of Bi(2212) show a broad stretching mode of B\ g symmetry at 
480 cm" 1 corresponding mainly to the vibrations of the apical oxygen O4 in the Cu-O-Bi 
bonds. The dc resistivity measurements of Crommie and Zettl [16] on Bi(2212) have been 
used to obtain values of e p and n from eq. (23). The sample with linear temperature 
dependence of the resisitivity gives n ~ 0.25 charge carriers per Sr atom and e p = 8 meV. 
For the Drude contribution we obtain [5] 

cj2 D = 47rn ef f e 2 /m" = 47rcr(O, T), (31) 

eff = ~^ nc e c h sech 2 (e p /?/2) (3 1 a) 

and 

r^ 1 =a 2 uj (3m*. (31b) 

The resistivity of Bi(2212) at 100 K is 120 u^-cm. This gives u pD = 6190cm" 1 and 
T D = 67 cm" 1 assuming m* nearly equal to the electron mass. These values are close to 
those obtained in ref. [14] from optical reflectance data. 

The expression for <?&(&, T] in eq. (30) applies for the correlated localized polarons 
when uj > 2e p . These act as oscillators having center frequency o>o which for Bi(2212) is 
480cm" 1 . There are localized as well as free polarons as shown in figure 1. The optical 
absorption for free polaron is given by eq. (15) while that for localized polaron is given by 
eqs (17) and (18). Therefore there are three types of contributions to optical conductivity; 
(i) Drude eq. (26), (ii) zero phonon scattering from localized polaron eq. (30) and (iii) one 
phonon (and higher phonon) scattering processes from free polaron states (eq. 15). 

The threshold for low frequency process (ii) in Bi(2212) is at 16meV which is 
observed by Reedyk et al [14]. The broad peak at 412cm" 1 arises both from (ii) and (iii) 
and the one at 996 cm" 1 arises from (iii) for a two-phonon scattering process which has a 
threshold at 960cm" 1 and a maximum at 996cm" 1 . 

4.2 Photoemission spectroscopy 

The angle resolved photoemission spectra (ARPES) of high- T c superconductors has not 
understood rlesrrite SftVftral attemnts R.7 81. Thp.rp. am however some, svstp.maHns 



that are emerging from these studies. In underdoped samples of I^SriCaCi^Og+tf a 
broad feature is always present at 100-200 MeV in the spectra taken near (TT, 0). Recently 
Shen and Schrieffer [3] have found that this feature is absent for overdoped samples along 
FS and XS directions of the Fermi surface (FS) while it is absent only along FS direction 
in the underdoped samples. They believe that this is due to strong dressing of the 
photohole for k ~ (TT, 0) in the underdoped samples. 

The difficulty in the interpretation of the angle-resolved photoemission spectra of high 
T c superconductors has been discussed by Mott [17]. Generally ARPES enables a Fermi 
surface to be mapped and the results are similar to band calculations. However in a doped 
semiconductor where holes near the Fermi surface drop down and form bosons there are 
no states available near the FS, so no excitation occurs and the spectrum does not map the 
FS. Following Mott we then look for the dynamical feature of the charge carriers in the 
ARPES spectrum of the high T c superconductors. 

The localized hole states are below the bottom of the polaron band by 2e p as shown in- 
figure 1. This is 16meV in Bi (2212). The bottom of the polaron band is below the atomic 
level by HUQ 60 meV in this system. In the photoemission spectroscopy the one-phonon 
optical absorption band following eq. (15) should peak at u> = UIQ = 60meV. But as in the 
ground state the polarons are in the localized state below the polaron band, the one phonon 
line appears at 16 meV below the Fermi surface. The line shape for the quasi-like excitation 
is given by eq. (15). In figures 4(a,b, c) we have plotted the intensity of absorption from 
eq. (15) with UQ = 60meVand a = 2.5, 1, 0.75 and compared them with the ARPES data 
of ref. [3] for overdoped F-S, X-S and underdoped F-S directed free polaron energy 
distribution curves, respectively. It follows that the electron-phonon coupling is stronger in 
the overdoped compared to the underdoped sample. The polaron along X-S in the 
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Figure 4(a). Energy distribution curve along F-S symmetry direction in the 
Brillouin zone of Bi(2212). The experimental data is taken from ref. [3] for the 
overdoped samples. The theoretical curve is from eq. (15) with m*/m = (1 - a/6)" 1 
and Q = 2.5. The binding enerev is in terms of ontical nhnnnn p.nf.rav Zv,^ whm-p M ic 
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Figure 4(b). Energy distribution curve along X-S symmetry direction in the 
Brillouin zone of Bi(2212). The experimental data is taken from ref. [3] for the 
overdoped samples. The theoretical curve is from eq. (15) with m* /m = (1 a/6)" 1 
and a = I . The binding energy is in terms of optical phonon energy TIUJQ, where UQ is 
the frequency of Bi 3+ -04 Cu 2+ bond stretching vibrational mode of Bi g symmetry 
and is taken as 480cm" 1 from ref. [15]. 
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Figure 4(c). Energy distribution curve along F-S symmetry direction in the 
Brillouin zone of Bi(2212). The experimental data is taken from ref. [3] for the 
underdoped samples. The theoretical curve is from eq. (15) with m* /m = (1 ct/6)~ 
and a 0.75. The binding energy is in terms of optical phonon energy HUQ, where U>Q 
is the frequency of Bi 3 +-O4-Cu 2+ bond stretching vibrational mode of B\ g symmetry 
and is taken as 480cm" 1 from ref. [15]. 



underdoped is localized and its line shape is given by the contribution from eq. (30) with 
additional contributions from multi-phonon processes adding at binding energies, 



a Vo? , T k(" - M)), (32) 



where ($) is a step function 

= for < 



This accounts for the broadness of this line. Later we show that evidence of multi 
phonon processes exist in femtosecond absorption studies. 

In several studies of photoemission spectrum hole pockets at (TT, TT) have been reported 
which are consistent with the transport and optical data. For underdoped samples, 
Marshall el al [7] found in Bi(2212), clear Fermi level crossings along the F-S line but 
over small portions of the BZ and its Fermi surface does not extend to the BZ edge. 

In the present model with 04 in the Sr-O plane in Bi(2212) the hole pockets are centered 
at the Sr atoms. From figure 2 it can be seen that the charge transport from one Sr atom to 
a neighbouring Sr atom can take place along the XS direction and to O atom along the F-S 
direction if O and Sr are both in mixed valence states and charge transport is produced due 
to the modulation of the Modelung energy by optical phonons [11]. In the underdoped 
samples with less number of holes the polaron state, that move along the XS direction from 
one Sr atom to another, cease to exist in the normal state and the polaron along this direc- 
tion tend to localize. On entering the superconducting phase the polarons along this direc- 
tion become delocalized and connect the hole pockets centered at (TT, TT) in the unit cells. 

In figure 2 linear chains of (2~ 2 + 2~) and (1~ 1 + 1~..) ions are seen in (1 10) 
direction. These alternate chains exchange their charge states when coupled with phonon 
at a frequency comparable to the phonon frequency. As argued by Shen and Schrieffer [3] 
if a photohole is created in these chains it excites a collective mode along (TT, TT) to attain 
the equilibrium state. 

4.3 Femtosecond optical absorption 

Femtosecond pump-probe measurements performed on Y-123, Bi(2223) and non- 
superconducting, YBa2Cu3O6+ y (y < 0.4) for sample temperature ranging from 7 K to 
300 K have been reported by Chwalek et al [8]. The observations are made in 
transmission on thin films with a laser that produced pulses of ~ 100 fs duration at 100 
MHz repetition rate and a wavelength of 630 run (1.98 eV). The change in refractive 
index was measured by passing a second beam polarized perpendicular to the pump 
beam. The change in transmissivity was recorded versus the pump-probe delay time. 

The experiment relates to the observation of the cross-modulation produced by the 
thermomodulation process. We take that the modulation is produced by the polarons 
whose frequency, CUQ, of the order of 50meV is small compared to the laser frequency fi 
(~1 eV). For large fi, the refractive index n is given by 

u 2 

= 1 7, (33) 

(fi + w) 2 v ' 
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Figure 5. Transmissivity change normalized to unity is plotted relative to the pump- 
probe delay time for the Bi(2223) film at 8 K and 202 K. The experimental data is 
taken from ref. [81. The theoretical curves are obtained from eq. (15) with 
m*/m= (1 a/6)~ and assuming that (LJ/UJQ 1) oc the normalized delay time 
(t/to 1) (eq. (34)). a for the theoretical curve for the superconducting phase 
(T = 8 K) is 1, and nonsuperconducting phase (T = 202 K) is 4.5. 

where OJ P is the plasma frequency. It can be readily shown that 



-1 



t-tp 

to 



(34) 



where (t to} is the pulse delay time and to is the reference time taken as the time of 
maximum transmissivity. 

In figure 5 we show the transmissivity change normalized to unity plotted relative to the 
pump-probe delay time for the Bi(2223) film at 100 K from ref. [8]. The theoretical plot is 
obtained on the assumptions that at low energies only free polarons are excited but when LJ 
exceeds 1.5 UJQ localized polarons are also excited. We then can write with = UJ/UJQ, 

r(0) = r(f/f ) = c/ieb[ fll ri(0) +a 2 r 2 (0)], (35) 

where FI (9} is given by eq. (15) and ^(0) is given by eq. (17) with cr(w, T) given by eq. 
(30). We take #2 as zero for 9 < 1.5 and gradually increasing for 6 > 1.5, 
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for e> 1.5 



for e< 1.5. 



(36) 



a\ simply normalizes the contribution from free polarons to its maximum value near 
0=1 and is 0.68. There appears to be a good agreement for a = 4.5 for the free 
polarons. 
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A similar analysis or me spectra ox me same sample at s K wnen it is superconducting 
shows that a\ =ai= 1 and a is 1.0 in this state. Consequently there is a delay in the rising 
edge near Ops and the line is much broader as shown in figure 5. Since, in eq. (16), 
(ft/2mu;o) is a measure of the size of the polaron, the decrease in a from 4.5 to 1 
indicates that in the superconducting phase the free polaron size increases by a factor of 
more than four compared to the normal state to enhance the overlap and increase the 
pairing interaction. 

The evidence of multiphonon like processes, as discussed in eq. (32) comes from the 
transmissivity change observed in ref. [8] for the nonsuperconducting YBCO sample at 
temperature ranging between 7 K and 296 K. In each case the main peak is followed by a 
large number of small oscillations in the tail of the relaxation whose amplitude decreases 
with increasing time as well as with increasing sample temperature. As discussed by 
Feynman et al [18] such oscillations are expected whenever external field supplies an 
energy quantum equal to the energy HUJQ of an optical phonon plus the energy of 
excitation to higher levels of the polaron. 



5. Conclusion 

We have examined the optical absorption, photoemission and femtosecond absorption 
spectrum of high temperature oxide superconductors and showed that these can be 
explained on the basis that the charge carriers are polaronic in nature and are primarily 
located in the Sr/Ba-O planes. 
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bstract. This paper reports conduction mechanism in a-SbaSes over a wide range of 
mperature (238 to 33 8 K) and frequency (5 Hz to 100 kHz). The d.c. conductivity measured as 
function of temperature shows semiconducting behaviour with activation energy AE = 0.42eV. 
lermally induced changes in the electrical and dielectric properties of a-Sb2Sea have been 
:amined. The a.c. conductivity in the material has been explained using modified CBH model. The 
md conduction and single polaron hopping is dominant above room temperature. However, in the 
wer temperature range the bipolaron hopping dominates. 
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Introduction 

uring the last few years amorphous semiconductors have been used in the manufacture 
solar cells, batteries and phototransistors as well as in some steps of technological 
ocessing of very large scale integration (VLSI) microelectronics circuits [1]. Thermally 
duced structural and electrical effects are known to be important in inducing the 
emory phenomenon in semiconducting chalcogenide glasses [2-6]. Alzewel et al [7, 8] 
ive studied the electrical conductivity of powdered antimony chalcogenides. Antimony 
loys have recieved great attention in the past due to semiconductivity of Sb 
mcentration. It was reported by Mott and Devis [9] that the effect of Sb in Se is 
fen greater than in Te in promoting crystallization. Recently, study has been made on 
in and thick films of the system BiioSbjcSego-^ (x = 35, 40 and 45) [10]. Crystalline 
DiSes has been prepared and studied by several workers [11, 12]. Antimony triselenide 
and gap ~1.2eV), [13] a p-type semiconductor having an orthorhombic D| structure 
ith unit cell dimensions a = 11.77 1, = 3.962 7 and c = 11.62 1 A [14] is used 
. the preparation of recording optical laser disc and as an insulating material in MIS 
)lar cells. The preparation and study of amorphous-Sb2Se3 can throw light on the 
Dssible applications of this material in that form. Measurements of electrical properties 
: a-SbiSes are reported in this paper over a wide temperature and frequency 
nge. Results indicate that electrical properties can be explained on the basis of 
idely accepted modified-correlated barrier hopping (m-CBH) model for chalcogenide 
.asses. 



2. Experimental details 

Ingots of glassy material Sb 2 Se 3 were prepared by melt quench method. The 99.99% pure 
elements were sealed in an evacuated (1 x 10~ 6 torr) quartz ampoule (12mm diameter 
and 160mm long). The ampoule was heated in a muffle furnace at the reaction 
temperature (6303C for 24 h). The composition <Sb:Se: 1 : 1) of the glasses was 
determined by EDAX and amorphous nature was verified by X-ray diffraction pattern 
obtained using X-ray diffractometer (Philips PW 1130/60). The pellets of 0.677cm 
diameter were prepared in the hydraulic press at a pressure of w 10 7 Pa. 

The measurements were carried out in a specially designed vacuum cell in the 
temperature range 233 to 313 K and frequencies from 6 Hz to 100 kHz. A sophisticated 
computer controlled ac-impedance system (Model 5206, EG & G, PARC USA) was used 
for the measurements [17]. Polished samples with a thickness 0.04-0.10 cm and area 
0.36cm 2 were put in contact with aquadag (a conducting emulsion) in a sandwich 
configuration. The sample was kept in vacuum in a copper sample holder ensuring no 
temperature gradient between the electrodes and the sample. The effect of temperature is 
studied using a closed cycle refrigerator (Model F-70, Julabo) to obtain different 
temperatures, which can maintain constant temperature within 1 K for all the 
measurements. Ohmic contacts were confirmed through linear I-V characterstics in the 
voltage range. All electrical measurements of real and imaginary components of 
impedance parameters (Z' and Z") and real and imaginary component of admittance 
parameters (Y r and Y") were made over a wide range of temperature (233 to 313 K) and 
frequency (6 Hz to 100 kHz). 



3. Results and discussion 

3.1 Dipolar behaviour o 
Figure 1 shows the plot of lno- dc versus 1000/T which was obtained using the expression 

(1) 



The plots in figure 1 are found to be linear over the temperature range studied. In the 
above expression C' = <JQ exp(7/&) and 7 is the temperature coefficient of the band gap. 
From the figure, the activation energy AE and C' for the sample are found to be 0.43 eV 
and 3.35 x 10~ 4 n -1 cm" 1 respectively. It has also been observed that AE is constant for 
different pellets (thickness = 0.04cm to 0.08cm) with a statistical error of 0.3%. 

The variation of log |Z| versus log / (figure 2) indicates that the impedance is sensitive 
to frequency at lower temperatures and gradually becomes independent of frequency with 
increase in temperature. This type of behaviour indicates dominance of band conduction 
at higher temperatures. 

The sensitivity of Sb2Se 3 to temperature is shown in figure 3. The figure gives the 
Nyquist plots (Z' versus Z" ) for Sb 2 Se 3 at different temperatures over the frequency 
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Figure 2. Frequency dependence of Iog|Z|(fi) over the entire temperature range 
studied. 



reduces with rise in temperature and regain their original size when temperature is 
increased. 

The Nyquist plots for Sb 2 Se3 form a perfect arc of a semicircle with its centre lying 
considerably below the abscissa (i.e., distribution parameter a 0). Figure 3 also gives 
the value of a calculated at different temperatures. In dielectric materials, the finite value 
of the distribution parameter a and a depressed arc are typical for a dipolar system 
involving multirelaxation processes [15-17]. But the dependence of a on temperature 
indicates that the mean position of the dipole is varying with variation in temperature. 
The dipolar nature of Sb2Se3 is further confirmed by investigating the variation of 
capacitance with temperature at different frequencies. It is clear from figure 4 that in the 
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Figure 3. Nyquist plots (Z' versus Z") at different temperatures. Variation of angle 
a (1 s)7r/2 with respect to temperature. 
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Figure 4. Temperature variation of capacitance at different frequencies. 

temperature. However the rate of change of capacitance (i.e., dC/dl 1 ) is higher for lower 
frequencies and decreases with increasing temperature. These results can be explained by 
the fact that thermally assisted hopping results in increasing the capacitance of the 
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R-C NETWORK MODEL 
Figure 6. R-C network model for lossy capacitance. 

material. In other words, in the lower temperature range, the dipoles remain frozen and 
attain rotational freedom when the temperature is increased. Thus the relaxation effects 
are confirmed by the rate of increase in capacitance with frequency. 

The inversion of Nyquist plots in the admittance plane are shown in figure 5. It is clear 
from the figure that the angle of inclination is less than rr/2. This type of behaviour has 
been explained on the basis of R-C network model (figure 6), such that total impedance is 
given by [18] 

Y(u) = B(iu} n + iuC a + Go 

= Go + B{cos(j7r/2) 4- i sin(s7r/2)~K + foC a . (2) 



10 



-5 



10 



-6 



-7 



10 



10 r- 



2.8 



bipolaron +bingie poiaron 
bipolaron only 
jingle poiaron only 



* 10 kHz 

A 31.5 kHz 

* dc 




3.2 



3.6 



4.0 



4.4 



1000/T (K~ 1 ) 



Figure 7. CBH fitting - the variation of v ac (ft~ l cm" 1 ) with respect to 1000/7 for 
experimental points and theoretical curves (m-CBH) at different frequencies. 



3.2 CBH model for a- 

The complex impedance plots show the dipolar multi-relaxation nature of a-Sb2Se3. This 
behaviour is explained by the R-C network model (figure 6). However, this does not 
further elucidate the nature of conduction mechanism in this sample. Therefore a.c. 
conductivity a ac (o;) at different temperatures is calculated. It is found that dispersion in 
the temperature dependence of cr ac (o;) curves increases at lower frequencies while it 
merges at higher frequencies (figure 7). 
From the expression 

the a.c. conductivity is obtained, where <7d c is given by (1) and its behaviour is shown in 
figure 1. The cr t0 tai(w) is measured directly from lock-in-analyser (2). The u ac (w) 
dominates over a& at lower temperatures, while <7d c dominates at higher temperatures 
(figure 7). 

The modified correlated barrier hopping (m-CBH) model explains the experimental 
results reported in this paper. It states that bipolaron hopping between the charged defects 
states D + and D~ is responsible for the a.c. conductivity in these semiconductors. Dipolar 
multi-relaxation behaviour of this sample is already indicated from the Nyquist plots 
described earlier. So the charged defect states are expected to be present in a-J 
According to CBH model a.c. conductivity is given by [19,20] 






(4) 



Electrical properties of a-antimony selenide 

where n is the number of polarons involved in the hopping process, NNp is proportional to 
the square of the concentration of the states and e' is the dielectric constant. R^ is the 
hopping distance for the conduction ur \ and is given as 

R^ 

where W m is the maximum barrier height and TQ is the characteristic relaxation time for 
the material. 

Figure 7 shows the variation of a.c. conductivity with temperature at different 
frequencies. It is evident from the figure that <r ac is very sensitive to temperature in the 
higher temperature regime. The low temperature a.c. conduction can be explained by 
considering bipolaron hopping between D + and D~ centers whereas the higher 
temperature behaviour is due to thermally activated single polaron hopping. At higher 
temperatures a number of thermally generated D centers are produced with a 
temperature dependent concentration [21] 

where NT is the concentration of D + or D~ centers at T K. The defect concentration 
factor NNp in (4) is replaced by 

NNp=Nj/2 (for bipolaron hopping) , 

NNp = Nj/2Qxp(U & {f/2kT} (for single polaron hopping). 

Thus for single bipolaron hopping this factor is thermally activated and hence cr ac (o;) is 
also activated. 

According to CBH model [22] 

B = W m - Wi + W 2 , (8) 

where B is the optical band gap, W m is the maximum band width and behaviour with 
frequency u is as 

a ac (u>}=Au s . (9) 

Here [19], 

dlnw 

while W = W m for bipolaron hopping and W = W\ or W 2 for single polaron hopping. 
Therefore, W m is estimated from the s values at lower temperatures. Values of Wi, W 2 and 
t/ eff were estimated from the values of s at higher temperatures. These parameters are 
then fed to CBH model to fit the a.c. conductivity data calculated from the experiments. 
The behaviour of (10) (s = 1 - 6kT/W*) can be seen in figure 8. This indicates that W* 
follows W m for lower temperatures and W 2 for higher temperatures. Figure 8 also shows a 
sudden change in conduction mechanism at a certain temperature (27 6-21 S K). Similar 
K^avinnr ie nntirpH in the. admittance dots Cfieure 5X The different parameters used for 
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Figure 8. Temperature dependence of j(= 1 - 6kT/W*) for theoretical curves (m- 
CBH) and experimental points. 
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Figure 9. Variation of./uA) with respect to temperature (K) at different 
frequencies. 

Further, the effect of temperature and frequency on the hopping length R U) which is a 
measure of the effective length of a dipole can be studied using CBH model. The 
variation of R u with frequency at different temperature is shown in figure 9. The figure 



indicates that /? w is more sensitive to frequency in the high temperature regime and the 
sensitivity decreases with decreasing temperature. Thus at low temperature, R^ is 
constant with frequency, which implies a constant value of capacitance which is found to 
be true from figure 4. 

4. Conclusion 

Nyquist plots or complex impedance studies confirm the dipolar nature of 802863 and 
multi-relaxation behaviour as seen in most chalcogenide glasses. These plots also indicate 
that the capacitive nature dominates over the resistive nature of the sample at lower 
temperature (<273 K ), while the resistive nature is dominant at higher temperatures. The 
conduction mechanism has a sudden change from bipolaron to single polaron hopping 
(276-278 K). The temperature and frequency dependence of a.c. conductivity is well 
explained by the m-CBH model. The contribution to a.c. conductivity from single polaron 
hopping is dominant at higher temperatures. 
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Abstract. A new type of correlated wave function for the normal helium is proposed in the form: 
\&(ri,r2) = )c m $ m with $ m = exp[ Q(TI + r 2 )]/(&ri 2 + a) m where a,b, a. are non-linear varia- 
tional parameters. Optimizing these parameters by the Monte-Carlo technique, an energy eigen- 
value of 2.903645(a.u.) is obtained with only three terms in the basis expansion alongwith 
satisfactory cusp condition, compared to the essentially exact non-relativistic value of 2.903724 
as given by Freund et al (Phys. Rev. A29, 980 (1984)). 
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1. Introduction 

Since the early days of quantum mechanics the study of electron-electron correlation in 
atomic and molecular physics has been in the lime light for its vital importance and 
significance in determining the behaviour of the bound electrons. A large number of 
important contributions [1-7] has been made on two-electron atomic systems during the 
last six decades. Of these, Pekeris [5] made an outstanding contribution in which he 
predicted essentially exact eigenvalues and other parameters for the ground and first few 
excited states of two-electron atomic systems using the Raleigh-Ritz upper-bound and 
Rubinow lower-bound methods. When the number of electrons become higher in an atom 
or a molecule, electron-electron correlation is taken into account by the many-body 
perturbation theory (MBPT), the configuration interaction (CI) method, the multi- 
configuration self-consistent field (MCSCF) theory, etc [8-10]. In all these methods, the 
calculations and the related computations become enormously difficult and complex. 

Recently some significant results have been reported in the literature on the two- 
electron atomic ground and higher excited states [11-13]. Umrigar et al [11] used Monte- 
Carlo multi-configuration self-consistent field theory with several hundred parameters in 
the basis set to consider electron-electron correlation in the normal H~, He, Be and few 
of their excited states. On systematically improving the wave function they were able to 
obtain accurate results in accord with Pekeris [5] and Freund et al [16]. Bhattacharyya 
et al [12], on the other hand, used a very simple three-parameter correlated wave function 



- Aexp (-/^ 12 JJ (i) 

and obtained an analytic expression for the energy eigenvalue of the normal helium: 
e= (t|#|#)/{#|#). They reported a value of -2.8894 (a.u.) (-78.6253 eV) which 
is much improved over the corresponding Hartree-Fock (HF) result, 2.86168(a.u.) 
(-77.87 104 eV), for the following parameters : a = 1.8395, A = 0.586, /z = 0.379. They 
obtained a value, 0.5365, for the correlation cusp condition, ((l/^)(d^f/dri2)) ri2 _ Q , as 
compared to the exact value, 0.5. Regarding an over estimation of the cusp condition by 
about 6%, it is possible that proper optimization of the parameters might reflect on its 
actual value corresponding to this wave function. Very recently, Drake [13] performed an 
elaborate calculation with the inclusion of non-relativistic and relativistic effects and used 
a huge basis set to report very accurate results correct up to 9 to 11 significant figures. 
In this paper , we propose a very simple form of the correlated wave- function for the 
ground state of the helium atom in the form 



c m * M (2) 

with 

$ m = exp[-a(n + ii)]/(fcr, 2 + a)", (3) 

where c, n 's are the linear expansion coefficients (m = 1,2, ... ,AT) and a,a,b are non- 
linear variational parameters. We use the Raleigh-Ritz variational method and exploit its 
upper-bound principle to obtain our results. The non-linear variational parameters are 
optimized by using the Monte-Carlo technique. 

The plan of the paper is as follows. We give a brief account of the theory in 2 and 
discuss the numerical method in 3. The results are discussed in 4. Atomic units (a.u.) 
are used in our calculations in which a = /no = e = ?i=l,l a.u. of energy (= e 2 / a o) = 
27.21 165 eV. 

2. Theory 

The standard eigenvalue equation according to the Raleigh-Ritz variational principle for 
the determination of the correlated orbitals $ m (ri,r2) takes the form 

(4) 



where the Hamiltom'an matrix H_ and the overlap matrix S have the elements H nm , S nm 
respectively over the functional space of the electron co-ordinates ri , r 2 (m, n = 1 , 
2, . . . , N}, the vector c' = (cj, C2, . . . , c) and e, the required energy eigenvalue. 

In order to evaluate the matrix elements H nm and S nm , use is made of the Fourier 
integral transformation for functions of the type exp( Ar)/r in the form 

1 f e ' P ' F rdp, (A>0) (5) 



and the Laplace integral transform for function (br\ 2 + a) m in the form: 
1 i r 00 ,. 



Correlated helium wave junction 



where F(m) denotes the gamma function, F(m + 1) = mF(m),F(l) = 1. Final expres- 
sions involve only a single-dimensional integral over [0, oo) which is evaluated 
numerically using 20 Gauss-Legendre quadrature points. Typically, the overlap matrix 
elements S nm are obtained as 



S = 



(7) 



with 



3. Numerical evaluation and Monte-Carlo optimization 

We have solved the eigenvalue equation (4) by first premultiplying it with an inverse of 
the matrix S and then finding the characteristic polynomial corresponding to the matrix 
S~ 1 H = A, say. The solutions of this polynomial constitute the required eigenvalue 
spectrum, the least of which gives the desired eigenenergy. For this we require a range of 
values for the non-linear variational parameters a, b, a, which is accomplished with the 
Monte-Carlo optimization technique. 



3.1 Krylov's algorithm 

In order to obtain the characteristic polynomial in e, we have used the Krylov algorithm 
for expanding the secular determinant, det(e/ A) [14]. Let 

det(e/-A) =e n +p 1 e"- 1 + +p n - l e+p n (8) 

be the characteristic polynomial. By using now the Cayley-Hamilton theorem which 
states that, the matrix A reduces its characteristic polynomial to zero, we obtain 



A" + piA n ~ l + + p n -iA+p n I = 0. 
Let us consider an arbitrary non-null vector: 



,() 



,(o) 



(9) 



(10) 



Postmultiplying both sides of equation (9) by y$\ we get 



,.(0) 



fii\ 



T ' -rPlT ~1 \-pn-\r ' -T AW M 

which takes the form 



Tr' J 2 " 5 




/^1\ 

P2 




3* 


n-l yi-2 ... 

\-'fl ^M '/I / 




U'J 




L? 1 / 



where 



(14) 



,(*) 



v??'/ 



(15) 



Since on the basis of the formula (12) yW = A/" 1 , (* = 1, 2, . . . , /i), the coordinates 
yf\y*2\ . . . ,yfi of the vector y_^ are successively computed from the relation 



(1) 



(0) 



(2) 

y - 



/=! 



(n) 



,,( M ~ 1 ) ; - 



(16) 



These are the input elements for the coefficient matrix in the system of linear 
equations (14) for the determination of the coefficients pi,p2, ->Pn of th e poly- 
nomial (8). It is to be noted that the coordinates of the initial vector y_W of eq. (10) are 
arbitrary. If the system of equations has a unique solution, then the roots pi,pi) ,Pn 
are the coefficients of the characteristic polynomial (8). If the system (14) does not have 
a unique solution, the problem gets complicated and a new initial vector is to be tried. In 
our calculations, we have used the unit vector (1, 0, . . . , 0) r as the starting solution. We 
have solved the linear system of equations (14) by using the scheme of Khaletsky [14]. 
The roots of the characteristic polynomial is next determined by using the Newton- 
Raphson method for a set of values of the non-linear variational parameters a, b, a. 



3.2 Monte-Carlo optimization 

Optimization of the non-linear variational parameters by the Monte-Carlo technique is 
another noteworthy point of our calculation which helps determine the energy eigenvalue 
with greater ease. Though use of this technique may have been made in bound-state 



state problems. We generate random numbers within a prescribed range of values for each 
of these variational parameters. Our philosophy of optimization of the non-linear 
variational parameters has been a departure from the conventional approach. We do not 
increase the number of basis terms in the wave function for a pre-determined set of values 
of these non-linear variational parameters. Instead, we first decide on the number of terms 
that are to be retained in a particular calculation. This is done with a view to fully utilize 
the potential of the wave function with the help of the Monte-Carlo technique. The non- 
linear variational parameters a,b,a are changed in a random manner within their 
prescribed ranges of variation one after the other, until the best energy is obtained by 
utilizing the upper-bound principle of the Raleigh-Ritz method. 

In order to check this procedure, we have also studied the following wave functions: 



, r 2 ) = ci exp[-(An + Bn)} + c 2 (i <-* 2) , (17) 

,r 2 ) = #exp[-o:(ri + r 2 )}[ Cl + c 2 exp(-/zr 12 )]. (18) 

The wave function (17) is known after Eckart [15] and incorporates radial correlation 
through the exchange symmetry. Results of the optimized parameters and the eigenenergy 
are listed below alongwith other values available in the literature. 

4. Results and discussion 

The present results are listed in tables 1, 2 and 3, alongwith the accurate values of 
Pekeris [5], and Roothaan and Weiss [6]. It is a noteworthy feature of our wave function 
that reliable results are obtained with only three terms in the basis expansion (2). The 
reason for such a finding is believed to be the inclusion of electron-electron correlation 
through higher orders by considering terms of the form l/(br\2 + d) m in the basis wave 
function. 

Furthermore, the basis functions <& m (r\ , r 2 ) of eq. (3) remain finite at the origin as r\, 
YI, r\2 all tend to zero. As is further desired, these correlated orbitals also vanish far away 
from the origin, when both r l5 r 2 >oo, r 12 remaining finite or becoming very large. As 
the bound-systems consisting of the atomic or molecular particles is confined in a finite 
volume in space, there arises no possibility of allowing r\i >co without making 
ri , r 2 > oo at the same time, which would any way make $ m (r i , r 2 ) > in the limit. Thus 
our choice of the correlated basis functions (3) seems to be rested on meaningful physical 
consideration. 

As mentioned earlier, we have also studied the same problem with different forms of 
the wave function as given by eqs (17) and (18). The Eckart wave function (17) yields a 
better value -2.875661(a.u.) for helium atom than the HF result -2.86168 (a.u.), while 
the present wave function (18) gives an eigenenergy -2.891 124(a.u.). As mentioned 
earlier, the Eckart wave function (17) represents a fair amount of correlation due to 
exchange symmetry. It is apparent from the values of the screening constants A and B that 
while one electron is assumed to be close to the nucleus, the other electron remains far 
away from it. Such a wave function which incorporates an additional flexibility for the 
electrons to avoid each other is often called a split-shell wave function. 



Table 1. Eigenenergies (a.u.) 
sion coefficients for correlated 



, the non-linear vanational parameters and the expan- 
wave functions alongwith cusp conditions. 



Form of the wave 
function (Z = 2) 



(a.u.) 



Parameters 



= 0.5 (Exact) 



1. Wave function (17): 

(a) Our work 

(b) Pauling and 
Wilson (1935) 

2. Wave function (18): 
(a) Our work 



(b) Roothan and 
Weiss (1960) 

(c) Bhattacharyya 
et al (1996) 

3. Wavefunction (2): 
(a) m = 2 



(b) m = 3 



-2.875661 A- 1.188, 5 = 2.184 

d = 1.0026, c 2 = 2.0053 
-2.8754 A = 2.15, B = 1.19 



-2.891124 a = 1.848, p = 0.01902, 0.372 

c, = 1.0, c 2 = -0.9513, 

N = 89.388 
-2.89115 a =1.8395 0.379 

-2.8894 a = 1.8395, p = 0.379 0.5356 

Cl = 1.0, c 2 = -0.586 

-2.891566 a = 1.855, a = 1.255, 0.376 

= 0.0001, c, =9644.931, 

Q = -12101.823 

-2.903645 a = 1.845, a = 1.935, 0.475 

b = 0.01, ci = 150.838, 
c 2 = - 132.873, 
c 3 = -298.482 



4. 'Exact' results: 


(see text) 






(a) Pekeris (1962) 


-2.903724375 


1078 terms 


_ 


(b) Roothan and 


-2.90039 


64 terms 


0.465 


Weiss (1960) 


(closed shell) 






(c) Freund et al (1984) 


-2.903724 


- 


- 



That the optimization procedure is numerically stable is demonstrated in tables 2 and 3, 
where we enlist the values of the eigenenergy for a number of non-linear variational 
parameters near the respective best value corresponding to the wave functions (17), (18) 
and (2) alongwith the averages for the nuclear and Coulomb interactions (l/r\} t (l/ni) 
correct upto sixth decimal place. The computer code is freely available for execution in 
personal computers using Fortran 77. The steadiness of the optimized values clearly 
emphasizes the correctness of our Monte-Carlo optimization procedure used in the 
evaluation. 

Further indication of the flexibility of our wave function is obtained by the values of 
the correlation cusp condition ((l/^)(d'^/dri2}) ri2 __ >0 . Compared to the exact value of 
0.5, we have obtained a value, 0.475, which is comparable to the result, 0.465, as reported 
by Roothaan and Weiss [6]. That the exact result is underestimated by the present wave 
function is quite natural. As the eigenenergy is symmetrically improved towards the exact 
value in the non-relativistic limit, the cusp condition is expected to gradually tend to its 
exact value of 0.5. This is evident from the figures of table 1 for the cusp condition 



and 



+ c 2 exp [- 



(18) 



Wave function (17): 
A B 



-e(a.u.) -e(eV) 



1.184 


2.184 


2.875654 


78.251293 


1.68400 


0.990102 


2.001885 


1.186 


2.184 


2.875660 


78.251456 


1.68500 


0.991048 


2.000824 


1.187 


2.184 


2.875661 


78.251483 


1.68550 


0.991521 


2.000293 


1.188 


2.184 


2.875661 


78.251483 


1.68600 


0.991994 


1.999762 


1.189 


2.184 


2.875660 


78.251456 


1.68650 


0.992466 


1.999231 


1.190 


2.184 


2.875659 


78.251429 


1.68700 


0.992938 


1.998700 


1.192 


2.184 


2.875656 


78.251429 


1.68800 


0.993880 


1.997638 


1.188 


2.181 


2.875659 


78.251429 


1.68450 


0.991401 


2.001644 


1.188 


2.182 


2.875660 


78.251456 


1.68500 


0.991599 


2.001016 


1.188 


2.183 


2.875661 


78.251483 


1.68550 


0.991796 


2.000389 


1.188 


2.185 


2.875660 


78.251456 


1.68650 


0.992191 


1.999135 


1.188 


2.186 


2.875659 


78.251429 


1.68700 


0.992388 


1.998509 


1.188 


2.188 


2.875654 


78.251293 


1.68800 


0.992782 


1.997257 


1.188 


2.190 


2.875646 


78.251075 


1.68900 


0.993176 


1.996006 



Wave function (18): c\ = 1.0 
a fj, -c 2 



-e(a.u.) -e(eV) 



1.846 


0.01902 


0.9513 


2.891121 


78.672175 


1.687294 


0.973528 


2.002286 


1.847 


0.01902 


0.9513 


2.891123 


78.672230 


1.688265 


0.974118 


2.001144 


1.848 


0.01902 


0.9513 


2.891124 


78.672257 


1.689236 


0.974708 


2.000004 


1.849 


0.01902 


0.9513 


2.891122 


78.672202 


1.690207 


0.975298 


1.998865 


1.848 


0.01800 


0.9513 


2.891021 


78.669454 


1.694998 


0.981041 


1.994185 


1.848 


0.01900 


0.9513 


2.891123 


78.672230 


1.689347 


0.974830 


1.999892 


1.848 


0.01901 


0.9513 


2.891124 


78.672257 


1.689292 


0.974769 


1.999948 


1.848 


0.01903 


0.9513 


2.891124 


78.672257 


1.689181 


0.974647 


2.000059 


1.848 


0.01904 


0.9513 


2.891124 


78.672257 


1.689125 


0.974585 


2.000115 


1.848 


0.01905 


0.9513 


2.891124 


78.672257 


1.689070 


0.974526 


2.000171 


1.848 


0.01906 


0.9513 


2.891124 


78.672257 


1.689015 


0.974465 


2.000226 


1.848 


0.01910 


0.9513 


2.891124 


78.672257 


1.688794 


0.974222 


2.000448 


1.848 


0.01920 


0.9513 


2.891123 


78.672230 


1.688242 


0.973617 


2.001002 


1.848 


0.01930 


0.9513 


2.891120 


78.672148 


1.687693 


0.973014 


2.001553 


1.848 


0.01902 


0.9500 


2.891094 


78.671440 


1.692237 


0.977991 


1.996976 


1.848 


0.01902 


0.9510 


2.891122 


78.672202 


1.689938 


0.975475 


1.999297 


1.848 


0.01902 


0.9512 


2.891123 


78.672230 


1.689470 


0.974964 


1.999768 


1.848 


0.01902 


0.9514 


2.891124 


78.672257 


1.689001 


' 0.974451 


2.000240 


1.848 


0.01902 


0.9515 


2.891124 


78.672257 


1.688765 


0.974193 


2.000478 


1.848 


0.01902 


0.9520 


2.891120 


78.672148 


1.687577 


0.972894 


2.001672 



tine Virial theorem states that -{P.E.)/(K.E.) = 2.0. 



parameters do not seem to be optimized properly, since the eigenenergy as reported by 
them is much worse than that quoted in Roothaan and Weiss [6] and that obtained by us 
independently with our procedure in the present calculation. 

In summary, the results seem to indicate that a satisfactory calculation can be made 
with mrre.latftd wave functions of the form (31 for studies on a two-electron svstem. This 
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f The Virial theorem states that -(P.E.)/{K.E.) = 2.0. 



prescription may be an inexpensive alternative to the well-tested correlated wave func- 
tions containing Hyileraas coordinates which have produced amazingly accurate results 
over the years ([2-6, 11]). 



5. Conclusion 

We have reported our findings with a new type of correlated wave function for the 
ground state helium atom. The results indicate that this wave function is sufficiently 
flexible to take into account the effect of electron-electron correlation in atomic systems. 
This application for a small system for which definitive results are available from 
the works of Pekeris [5], Drake and co-workers [13] and Freund et al [16] shows 
the simplicity of our approach. Now, neither the use of the form of correlation with 
inverse powers of r\i in the wave function such as eq. (3) is certainly new (see, for 
instance, Kinoshita [4] or Umrigar et al [11]) nor the optimization of the non-linear 
parameters by Monte-Carlo technique is novel. What we show in this work is that 
reasonably accurate values of the eigenenergy of the two-electron systems can hopefully 
be obtained with our prescription with ease and elegance. Such a simple form of the wave 
function may therefore find applications in other areas, such as, atomic or molecular 
collisions where bulk of the theoretical efforts are involved with inaccurate bound-state 
wave functions without taking account of any correlation of the atomic/molecular 
electrons. 



Further, this simple approach can be extended to write molecular wave functions 
through linear combination of atomic orbitals (LCAO) which would explicitly 
incorporate electron-electron correlation. Thus this method might be used alongwith 
other methods in structure calculations, such as, the density functional theory ([17-21]), 
the CI method, the MCSCF theory, the MBPT, etc. 

All these calculations were performed on an IBM-compatible PC/AT of 100 MHz 
clockspeed with double precision arithmetic. Further applications of this correlated basis 
function would hopefully make a new beginning in the study of many-body systems in 
atomic, molecular, nuclear and solid-state physics. 
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The scattering of low energy electrons by atoms and molecules has been studied for many 
years. It has been found to be important for the better understanding of weakly ionized 
plasma. The electron collision studies of nitrous oxide molecules are of particular 
importance as N2O lasers are now being used as a secondary frequency standard in areas 
of spectroscopy [1]. Also, nitrous oxide plays an important role in stratospheric chemistry 
and eventually affects the Earth's climate [2]. 

Total electron scattering cross sections for nitrous oxide were measured many years 
ago by Briiche [3] at electron energies larger than 1 .0 eV and by Ramsauer and Kollath [4] 
at energies varying from 0.1-1.2 eV. Since then, only three more measurements have been 
reported using high energy resolution electron spectrometers. Zecca et al [5] used a 
retarding potential difference electron monochromator to obtain normalized total cross 
sections at electron energies ranging from 0.75-7 eV whereas a linear transmission 
method was employed by Szmytkowski et al [6] to obtain absolute cross section values in 
the energy range 0.8-40 eV. Absolute total cross sections for 1 -500 eV electron scattering 
of nitrous oxide were also reported by Kwan et al [7] by a beam transmission method. 
However, some relative transmission measurements have been carried out by Boness et al 
[8] and Sanche and Schulz [9] in the energy range 0-6 eV, but unfortunately, it is not 
possible to obtain cross sections from these results. Also, no theoretical computations 
have so far been reported in the literature for obtaining the absolute total electron 
scattering cross sections for nitrous oxide. In a nutshell, there are only four measure- 
ments available for a possible comparison of cross section values at electron energies 
from 0.7-10 eV. Although the general shape of the cross section curves reported by 
Szmytkowski et al [6], Kwan et al [7], Zecca et al [5] and Briiche [3] is more or less 
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similar with respect to electron energy, the cross section values in some cases are 
different. There is a large discrepancy in cross section values at and around the well- 
known 2 ^2 + shape resonance centered at 2.3 eV energy [10]. Around the peak of the 
shape resonance, the cross section value reported by Szmytkowski et al [6] is about 26% 
higher than that obtained by Briiche [3] and is larger by a factor of 2.3 than the norma- 
lized cross section given by Zecca [5]. However, there is satisfactory agreement in cross 
section values reported by Kwan et al [1] and Szmytkowski et al [6] at the peak of the 
shape resonance. Also, there is a large discrepancy in the results obtained by the four 
groups on both sides of the shape resonance at energies ranging from 1.5-3.5 eV. In view 
of this, we decided to measure the absolute total scattering cross sections for nitrous 
oxide at low electron energies. 

This paper is part of the ongoing research programme taken up to measure total 
electron scattering cross sections for atoms and molecules at low electron energies using 
a photoelectron source. Previously, absolute total cross sections for helium, neon [11], 
argon, krypton, xenon [12], molecular hydrogen [13] and molecular oxygen [14] have 
been measured at electron energies from 0-1 OeV. In this paper, we present measurements 
of scattering cross sections for electrons scattered by nitrous oxide with projectile energy 
varying from 0.73-9.1 4 eV. 

A detailed description of the experimental set-up used for the study, the method 
for analysis of the data and error analysis has been reported previously [11,12]. In 
brief, the experiment consists of measuring the intensities of the peaks in the photo- 
electron spectra of the source gas such as argon, krypton, xenon or nitrous oxide itself. 
The photoelectrons thus produced by the source gas are scattered by the target gas, the 
cross section of which is to be determined. Each peak in the photoelectron spectrum 
of the source gas provides one point in the electron energy scale. Complete scanning 
of the electron energy is performed by varying the energy of the ionizing radiation as 
well as changing the source gas. Using different combinations of photons of three 
different wavelengths (Hel, 58.4 nm, Nel, 73.6nm and 74.4 nm) and four different 
source gases and neglecting the energy points where the statistics for the photoelectron 
intensity were poor, it was possible to measure scattering cross sections at 19 electron 
energies. 

The electron scattering cross sections for nitrous oxide were measured using the 
method described previously. When source and target species are different, the electron 
scattering cross section a can be evaluated by using the following equation: 



Here I e \ and / e2 are the amplitudes of the photoelectron peaks at two different gas 
pressures PI and P 2 , /AOI and I\ Q2 are the incident photon intensities monitored during the 
experiment by the beam monitor at the two pressures, l\ is the distance from the centre of 
the beam splitter to the circular aperture covering the ionization region, / 2 is the distance 
from the same aperture to the actual ionization region defined geometrically by different 
slits in the acceleration region and cylindrical mirror analyser, x is the scattering path 
length, a is the ratio of the pressures of the target gas outside and inside the ionization 
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Figure 1. Total electron scattering cross section for nitrous oxide as a function of 
incident election energy from 0.73-9. 14 eV obtained by different researchers. The 

symbols denote as follows: V, Briiche [3]; , Zecca [5]; O, Szmytkowski et al 

[6]; D. Kwan et al [7]; *, present work; , cross section curve after proper 

curve fitting (present work). 



have been taken from Samson and Yin [15]. The ratios of /^oi and /AOI were determined 
from the beam splitter as described previously by Kumar et al [11]. 

When source and target gas species are the same, the electron scattering cross sections 
could be calculated by the equation given below: 



Ini^ 

^e\ 



Pi - 



760 



(2) 



The electron scattering cross sections for nitrous oxide were measured using equations 
(1) and (2) respectively in the two cases described above. In both these equations, all 
parameters I e \, l e i, /^oi, I\Q2-> PI and PI could be determined experimentally and cross 
sections could be deduced. 



All errors in the measurement of electron-scattering cross sections have been discussed 
in detail previously by the present authors. In the present experiment, the most probable 
error was estimated to be 3%. The problem of inadequate discrimination against 
forward scattered electrons has been attempted to estimate the magnitude of this effect on 
the measured cross sections in the experimental set-up used by us. A Monte Carlo method 
was developed and the electron trajectories originating from the ionization region and 
terminating at the detector were numerically studied using five randqm generators i.e. the 
finite length of the ionization region where photoelectrons are produced, random angle of 
projection within the allowed cone of ejection determined by the geometry, spread in 
energy of the photoelectrons determined by thermal broadening, place where a chosen 
fraction of the emitted photoelectrons are undergoing forward scattering and deflection in 
angle within a specified cone due to forward scattering. Assuming a differential scattering 
cross section of 100 A 2 per steradian, the error in total cross section due to lack of 
discrimination of forward scattered electrons was found to be dependent on the overall 
cone angle of emergence. For our experimental set-up where the overall cone angle 
of emergence has been estimated to be 6, the maximum error has been found to be 
about 0.5%. 

The total electron-scattering cross sections for nitrous oxide as measured in the present 
experiment are shown in figure 1 for electron energies ranging from 0.73-9. 14 eV along 
with error bars at three energies only. Also shown in the figure are the measured cross 
sections of Briiche [3], Zecca et al [5], Szmytkowski et al [6] and Kwan et al [7]. The 
cross section values reported by Ramsauer and Kollath [4] is not shown in the figure as 
there are only 10 cross section values available at energies from 0.1-1.2eV and 3 values 
between 0.7 to 1.2eV. The energy width for electron beams used by Kwan et al [7] was 
about 100 meV and the experimental error for absolute total cross section measurement 
was about 5% at 1.2eV and 3% at higher electron energies. 

Szmytkowski et al [6] reported an overall experimental error (calculated as the linear 
sum of all single contributions) of 10% at the lowest electron energy, rising to 12% at the 
sharp maxima and decreasing to 8% at and above lOeV. The energy spread of the 
electron beam used by their group was about 70 meV. Relative total scattering cross 
sections were measured by Zecca et al [5] in an attempt to locate broad resonance, if any, 
and to determine the absolute size of the sharp resonances observed by Sanche and 
Schulz [9]. However, their relative cross sections were normalized at 4 eV to the value 
obtained by Briiche [3] at the same energy. Also shown in figure 1 is the cross section 
curve after proper curve fitting to the results obtained in the present experiment. This has 
been carried out using an appropriate mathematical function in the library of large 
number of such functions available in the Tablecurve-2D (Jandel Scientific) software for 
which the regression coefficient is closest to unity. It may be pointed out here that the 
theoretically computed electron scattering cross sections for nitrous oxide have not been 
shown in figure 1 as these have not been reported so far in literature by any research 
group. 

For the sake of discussion, the cross section curve (figure 1) in the entire electron 
energy range from 0.73-9. 14 eV has been divided into two regions. These include regions 
below 5.5 eV and between 5.5 and 9.14 eV. In the first energy region a sharp increase in 



involved was the 2 ]T f me N 2 ion. From figure 1, it is clear that the cross section 
values at the maximum of the resonance peak and also at the wings on both sides of the 
peak show large discrepancies. To have a clear comparison, the cross section curves 
reported by Briiche [3], Zecca et al [5], Szmytkowski et al [6] and Kwan et al [7] along 
with the curve obtained in the present experiment have been made to undergo proper 
curve fitting procedure. This has been done with the help of Tablecurve-2D (Jandel 
Scientific) software. The energy and the cross section value at the peak of the shape 
resonance thus obtained are 2.3 eV and 28.8 A 2 [6], 2.3 eV and 29.1 A 2 [7], 2.1 eV and 
12.5 A 2 [5], 2.3 eV and 21.9 A 2 [3], and 2.2 eV and 23.2 A 2 (present work) respectively. 
The electron energy value at the peak of the shape resonance has been reported by 
different research groups to be the same but the cross section values have been found to 
be different. The cross section values reported at the peak by Szmytkowski et al [6] and 
Kwan et al [7] are same within the stipulated experimental errors but these values are 
larger by about 25% as compared to the cross section obtained in the present measure- 
ment. The cross section value measured in the present work is about 6% higher than that 
obtained by Briiche [3] whereas the value given by Zecca et al [5] is about 46% lower 
than that reported in the present work. It has been mentioned earlier in the text that the 
relative data of Zecca et al [5] has been normalized at 4.0 eV to the cross section value 
reported by Briiche [3]. It is difficult to visualize why the normalization was carried out at 

4.0 eV because this energy neither corresponds to the maximum nor the minimum of the 
resonance peak. Also, the ratio of the cross section values at the maximum and the 
minimum (9. 1 A 2 at 4.4 eV electron energy) of the resonance peak after normalization is 
extremely small i.e. about 1.4. In comparison, such a ratio varies from 2.1 to 3.2, the 
smallest and the largest being obtained from data reported in the present experiment and 
by Kwan et al [7]. Incidentally, the shape resonance passes through the minima at 5.4eV, 

5.1 eV, S.OeVand 5.2eV having cross section values 9.6, 9.0, 8.3 and 11.0 A 2 as given by 
Szmytkowski et al [6], Kwan et al [7], Briiche [3] and the authors in the present work 
respectively. Except for Zecca's results, all other four measurements are fairly reliable. 
However, there are two sets of cross section values available from the four measurements 
for the peak of the shape resonance. One set of values, from Szmytkowski et al [6] and 
Kwan et al [7] is about 28% higher than that given by Briiche [3] and authors in the 
present work. At this stage, it may be worthwhile to discuss about error in all the four 
measurements associated with incomplete discrimination against projectiles at small 
angles in the forward direction. In the present experiment, the maximum error because of 
lack of such discrimination has been found to be about 0.5% whereas such an error in the 
measurements by Kwan et al [7] is also small, of the order of 0.5 to 1 %. Briiche [3] has 
not discussed about this potential source of error associated with their instrument, 
whereas Szmytkowski et al [6] have found an error of about 1 % at 60 eV and even smaller 
at lower energies. In view of the above discussion, it is clear that two sets of cross section 
values exist for the peak of the shape resonance and this controversy needs to be resolved 
in future. 

In the energy region from 5. 5-9.1 4 eV, the cross section curve after passing through a 
minimum on the decreasing wing of the main shape resonance peak shows an increasing 
trend in cross sections at higher energies. In the measurements reported by Zecca et al 
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case, the difference being more than the stipulated experimental errors. The cross section 
values reported in the present experiment in this electron energy region are the highest 
whereas those reported by Briiche [3] are the lowest. 

The values of scattering cross sections for nitrous oxide as measured in the present 
experiment are given in table 1 along with the photon wavelength, source gas and the 
photon state for the corresponding electron energy. 

As nitrous oxide and carbon dioxide molecules are isoelectronic and linear in their 
ground state, it may be interesting to find out the difference in the electron scattering 
cross section curves for the two molecules. The cross section measurements for the CO2 
have also been carried out by the present authors and these would be published elsewhere 
in detail. Our measurements in the case of these two molecules show two broad features. 
Qualitatively, the shape of the cross section curves of N2O and CO 2 has striking 
similarities, but quantitatively, it has been found that the cross section values for CO 2 are 
smaller than those obtained for N 2 O at all electron energies below lOeV. Secondly, the 
resonance peaks have been found at different electron energies ie. 2.3 and 3.8 eV for N 2 O 
and CO 2 respectively. It may be worthwhile to make an attempt to explain at least the first 
broad feature. N 2 is a polar molecule with a permanent dipole moment whereas CO 2 is 
a non-polar molecule. It has been found that at smaller angles, the differential cross 
sections at intermediate energies are larger in the case of polar molecules than those for 
the non-polar molecules [16,17]. It is presumed that the same is true at low electron 
energies also. In that case, the total cross sections would also be larger in case of 
molecules havine Dermanent dioole moment. 
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Abstract. The hyperspherical partial wave method (Das [7]) has been presented here in a 
symmetric form so that effects of coupling among different partial waves may be included in a 
systematic way. It is also outlined here how to solve the relevant coupled set of radial wave 
equations numerically. Some preliminary results are presented for S,P and D waves in the low 
energy domain of 30-50 eV for the incident electron. In this calculation only two important partial 
waves are included in each channel. The results are compared with experimental ones and appear 
very encouraging. 
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1. Introduction 

Three very successful methods of calculating total cross sections for ionization of 
hydrogen atoms by electrons at low energies are the convergent close-coupling (CCC) 
method (Bray and Stelbovics [1,2]), the /?-matrix with pseudostates (RMPS) method 
(Bartschat et al [3]) and the hyperspherical close-coupling (HSCC) method (Kato and 
Watanabe [4]). All these methods give good total cross values (see [4, 5]) down to very 
low energies. However their capabilities in reproducing the differential cross sections 
have not been tested (except in one or two cases such as in Bray et al [6]). Moreover each 
of these methods has its own difficulties and limitations. Another, very similar approach, 
is the hyperspherical partial wave (HSPW) method, suggested by the present author [7]. 
This is perhaps a better approach. Since this approach is straightforward and most natural, 
it is just an extension of the partial wave method so useful in two-body scattering 
calculations. It is also free from any pathological difficulty such as the non-existence of 
certain integrals for some of the angular momenta for the partial waves (see the second 
paragraph after eq. (2) of ref. [6]). It also does not require additional knowledge such as in 
HSCC calculation where one needs the values of the wave function in the asymptotic 
domain for matching purpose from some other sources. Moreover such values are known 
only approximately and no analytic solution exists for the asymptotic domain. Some 
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obtained following this approach, for Erhardt asymmetric geometry [8] and for coplanar 
constant 0\i geometry [9], 9\i being the angle between the outgoing electrons. In these 
calculations effects of coupling among different partial waves, in hyperspherical 
coordinates, are neglected. To include the effects of coupling among different partial 
waves it will be better to proceed in a symmetric form as is indicated below. In our earlier 
papers [7-9] we did not consider the problem of how to solve the resultant coupled set of 
radial wave equations numerically. Here we paid attention to this very problem and 
outlined a method. Following the present approach some calculations have also been 
done. Here we report also these preliminary results. In 3 we describe the present 
symmetric approach which intends to include fully the coupling effects. In 3 we 
describe the calculations of the present article. Results for S, P and D waves are reported 
in 4. Some concluding remarks are made in 5. 

2. Hyperspherical calculation in a symmetric form 

We consider ionization of a hydrogen atom in its ground state by an electron of incidence 
energy E, and momentum pi, the momentum of the ejected and scattered electrons being 
Pi and p2 and coordinates r\, r 2 . 

Then ionization amplitude for symmetry S is given by 



(i) 

where S = or 1 corresponding to a symmetric or an antisymmetric wave function. Here 




and ^f s ' is the exact scattering state wave function in the final channel with symmetry S 
and satisfying incoming converging boundary condition. Here we indicate how one can 
calculate ^ fs , and hence the scattering amplitude using hyperspherical coordinates 
(/?, Q,ri,r2), where 

o: = tan~ 1 (r2/ri). 




and use the notations u> = (a, n, r 2 ), cu = (ao,pi,p 2 ). 
Now we note that plane wave has the expansion [7] 

exp(i Pl TI+ *p 2 -r 2 )/(27r) 3 - (2/7r) 1/2 T i\j, x (p}/p 3/2 }^M^(^, (2) 



Henceforth A stands for the eigenvalue l\ +li+ 2n or the set (/j , / 2 , n, L. M) depending 
on the context and z/ A = A + 3/2, p = PR. For our present formulation it is necessary to 
consider the symmetrized plane wave function, which may be expanded in hyperspherical 
coordinates using the symmetrized angular wave function [10] 



(a)yr(n, *>), for /, = / 2 = 1. (3) 

These wave functions are orthonormal and have definite values of (L, M, S, TT), TT standing 
for parity and taking values or 1. Then one can easily prove that 

{exp(ip r ri + ip 2 -r 2 ) + (-l) s exp(p, -r 2 + zp 2 -ri)}/(27r) 3 

= 2(2/ 7 r) 1 / 2 V f A (7n(p)/P 3/2 )^(^o)^M. (4) 



We consider symmetrized final scattering state \l>[ s ' with incoming boundary condition. 
This state may be expanded as 

The radial wave functions F s x (p)'s satisfy the equations 

[(d 2 /dp 2 ) + 1 - (v x (v\ + l))/p 2 ]Fj + 2(27j A ,/p)Fj, = 0, (6) 

X' 

where 



C = (I/ cos a) + (I/ sin a) (l/|h cos a - r 2 sino;|). 

At this point an explanation is in order why the formulation of the problem in 
symmetric form is expedient. Because of conversation laws, for any fixed /^ = (L, 5, IT) 
only radial waves with this particular yu enter in a coupled set of equations, reducing 
considerably the computational problem. 

In the above equations those A, A' occur for which there are fixed values of L,M, S and 
TT. For the present case M = only needs consideration, since in the calculation of the 
scattering amplitude, given by equation (1), only states with M = gives a non-zero 
contribution. Henceforth we fix M = everywhere. Now we denote (L, 5, ?r) by the single 
letter and introduce a single index N in place of (l\, / 2 , n). Then we have Fj[(p) =f' J N (p} 
and equation (6) becomes 



[(d 2 /dp 2 )+ 1 -(KK + 1)/P 2 )]/P) + (27^'/P)A(P) = 0, (7) 



where v\ = 4 7^ A , = 7^, for each fixed p and N. In any realistic calculation N, N' are 
restricted to certain maximum value Nmx N I)VC depends on //, but this will not be expli- 
citly indicated here. 

Now we consider particular solutions of the above set of equations. It may be noted 
that there are two particular sets of solutions fm k \p) and t$ k \p) which behave 
asymptotically as 



00 
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and 

Here (3% is the fcth eigenvalue of the charge matrix A' 1 = (TJ^/) and 
fl 





are the fcth eigenvectors of it. We determine particular solutions f $ ' and f in the 
asymptotic region by substituting 



/=0 /=0 



where 0% = p + j3^ In 2p in the differential equation (7) and on equating coefficients of 
sin0jt/p /+1 and cos^/p /+l to zero. Thus we have the following recurrence relations for 
the unknowns 



-, (9a) 

- Z(2/ - l)/3?/]Yf ^ 
+ [(2Z - l)/3!jAf + /A - / 2 (/ - l)/]xf ~ l \ (9b) 

where 



The above equations with the initial conditions 

XM) _ vA 1 v^(^ n 
ft A * j x ft u 

give fjn '(p) and the conditions 



Now the physical solution will be a superposition 



k k 

The partial wave containing <p^ (u;) will be then 

(2/tr) 1 / 2 



AW corresponds to the set A for the index N and for fixed fj,. This must agree (apart from a 
Coulomb distortion) with the part 

2(2/7r) 1/2 * Aw exp(i(p - i4V 2 ))/ 2 ip 5/2 )^(wo)</^H 

of the symmetrized plane wave in order that *[~ ) be asymptotically a superposition of 
(distorted) symmetrized plane wave plus incoming converging waves. So we have, 

y^(c^ + id%)a% N = 2exp(i(A n - v N )ir/2)(p s (u ). (11) 

k 

Putting p' = c^ + id%,, q% = c^ id 1 ^ we have 

X^pM = _2exp(/7r/4) 
where 




Then 




(12) 



Here the inverse (X^)"" 1 exists since X^ is actually the similarity transformation which 
diagonalizes the symmetric charge matrix A^ . 

Next we construct N^ independent solutions of equation (7) which vanish at origin. 
We denote the corresponding solution vectors by f(> , k = 1,2, . .. ,N rnjc . Then the 
physical solution vector may also be written as 



We match values of radial waves (and their derivatives) given by right hand sides ol 



(14) 



where g RK and giK turns out in the forms 

gRK = E sS'rf'/z W + E &^N'I 



N 1 

and 



and v^(wo) and </?^, 7 (u;o) are the real and imaginary parts of ^(LJQ). In this way radial 
waves /Jv(p)'s are determined for the scattering states 



This on substituting in equation (1) gives the scattering amplitude and hence the cross 
sections. 

3. Calculations 

From the results of CCC and RMPS calculations of Bartschat and Bray [5] it is now 
known that at low energies S 1 , P and D waves give most significant contributions to the 
total cross section, P wave itself accounting for more than 50%. So we include in our 
calculations only the channel 1)3 5 e , li3 P and 1 ^D e . Here a channel corresponds to fixed 
values of (L, S } P} = /LI. We report two sets of results. In one of these only one, most 
important, partial wave in each channel is retained. For example 3 P includes the partial 
wave with J 1 =0,/2 = l,n = only. In the second calculation we include two important 
partial waves in each channel. These are shown in table 1. 

In our present scheme, determination of solution of the coupled set of equations [7] is 
the most difficult task. As is described in 2 the equations are solved by matching method. 
In the asymptotic domain solutions f f} (p) and f ^ (p) are determined from asymptotic 
expansion (equations (8) and (9)). Then another set of solutions are continued from the 
neighbourhood of origin, say starting at e (with e = 0.1). At this point N m solutions are 
constructed with initial values zero and derivatives SNK for the Mh element of the fcth 
solution vector. These are then continued numerically to pi nt , a point in the asymptotic 
domain. In this numerical integration of the equations (8) we adopted a Tailors expansion 
with 10 terms in steps Ap = 0.01 initially up to 0.2, then in steps Ap = 0.05 up to 0.2, 
then with Ap = 0.2 up to p,^, a point beyond which no contribution arises in the T-matrix 
elements. p- mt was varied from 20 a.u. to about 150 a.u. in steps of 5 or 10 a.u.. In general 
there was good agreement in the results for different matching points p int . But occa- 
sionally there were oscillations also. In some instances there were some spurious singular 
like behaviours. With further change in the value of p int , the results again become what it 



IQ 



3pO 



(0,0,0) (0,0,1) (0,1,0) (0,1,0) (1,1,0) (1,1,1) 

(0,0,2) (0,0,3) (0,1,1) (0,1,1) (1,1,2) (1,1,3) 



Table 2. Partial wave cross sections (in ofy for ionization of hydrogen atoms bj 
electrons in different channels in the energy range 30 50 eV for calculations witi 
A^c 2. Quantities in brackets represent those with N,, w = I . Total cross sections are 
also presented (which include cross section for F state of RMPS calculations). See tex 
for details. 
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0.06 


0.03 


0.07 


0.24 
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0.01 


0.53 


0.71 


0.69 



results calculated with different pj nt in the range noted above. Often average values wen 
taken as the final results. In these calculations scaling of the wave functions is necessar) 
since for a variation of p from say p = 1 .0 to p = 10 units the wave function changes b; 
many orders. As already stated we made calculations for N njx = 1 and A^ = 2 only. Fo 
further increase of N tnx often the results become unphysical. This implies that the 
stabilization of the calculations becomes essential. This problem is left for a futur 
investigation. Now in obtaining T-matrix elements all the integrations except the integrs 
tion over p were done analytically. The p-integration was done numerically. In calculatin 
the partial cross sections from the T-matrix elements the remaining integrations ar 
amenable to analytical evaluations. In this paper we present only the various partial cros 
sections and values of certain parameters related to these. 

4. Results 

Partial wave cross sections for ionization of hydrogen atoms by electrons in the ener 
range 30-50 eVof our present calculation are presented in table 2. In this table two sets < 
results are presented. One set corresponds to N nvc = 2 and for comparison a second set 
presented for A 7 ^ = 1. For 30, 35 and 40 eV energies the two sets of results are close 
each other indicating that at these energies most important contributions arises fro 
single states (0,0,0) in l S e , (0,0,1) in 3 S e , (0,1,0) in 1 P and 3 P, (1,1,0) in l D e and (1,1, 
in 3 D e . As the energy increases, other states in these channels become important. The tv 
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Figure 1. Spin asymmetry for electron impact ionization of hydrogen atoms as a 
function of the projectile energy. Theory: continuous, present calculation; dash-dotted, 
CCC [5]; dash-double dotted, RMPS [5]. Experiment: , Crowe et al [12]; o, Retcher 
etal [11]. 

states shown in table 1 in each channel give reasonable contributions to partial cross 
sections as well as to total cross sections in the energy range 30-40 eV. Total cross 
sections of the present calculation, the CCC and the RMPS calculations agree quite well 
(see table 2). Beyond 40 eV, however, the results go to depart considerably. We have not 
calculated F-channel contributions. F-channel gives about 5% to 10% contributions to 
the total cross sections at these energies. In calculating the total cross sections we have 
included the F-charmel values of RMPS calculation to our results. Next we look to table 3 
where we presented ratios of partial cross sections to the total cross sections. Our results 
agree well with the CCC results and with the RMPS results for 30 eV energy. For 50 eV 
energy our 5-wave results appear somewhat higher compared to those of other calcula- 
tions, while the D-wave cross sections are somewhat lower. The P-wave cross sections of 
these calculations which account for more than 50%, agree reasonably well with one 
another. 

Finally in figure 1 we have presented values of an important parameter of these cross 
section results. This is the spin asymmetry parameter defined by 

spin asymmetry = (a 5 - c?}!^ + 3a T }. 



et al [12]. It is interesting to note that our results for A^ = 2 agree beautifully with the 
experimental results and appear somewhat better compared to the CCC and RMPS results 
in the above energy range. Excepting one point of Crowe et al all these experimental 
points are very close to our theoretical curve. 

5. Conclusion 

The results of our present preliminary calculation of electron hydrogen atom ionizatior 
collisions, presented here in a symmetric form, and which includes fully the coupling 
effects, appear very encouraging. Thus simple inclusion of two most important states ir 
each channel gives good representation of the partial wave cross sections in the energy 
range 30-40 eV. However, for better results and for other energies, one has to include 
larger number of states in each channel. Matching method also does not appear verj 
efficient. Improved calculations including a large number of partial waves in eacl 
channel now appears relatively difficult. Then stabilization of the calculation for the 
solution of the coupled set of equations will be necessary. Once this is achieved man] 
interesting results and insights in the low energy domain are expected. 
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Abstract. A one dimensional numerical model to simulate the laser matter interaction in the 
presence of an ambient gas is presented. The model is developed by making appropriate 
modifications in MEDUSA, a one dimensional Lagrangian computer code, which simulates laser 
plasma interaction in vacuum. Various parameters of the plasma such as velocity, electron 
temperature, ion temperature, density, pressure, shock wave intensity of the plasma as it expands 
into a background gas are simulated. The results are compared with the experimental observations. 
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1. Introduction 

The study of laser-matter interaction has been a subject of paramount importance in 
recent years. There have been reports on laser-matter interaction in the presence of an 
ambient gas [1-6], electric field [7, 8] and magnetic field [9]. The genesis of this study is 
to understand the interaction of a high power pulsed laser beam with a solid target in the 
presence of an ambient gas and to develop a comprehensive numerical model for 
the process. The laser ablation at moderate intensity is particularly important because of 
its application in thin film deposition [3] and laser drilling [10-12]. Further, the 
interaction between laser produced plasma and ambient gas is of interest because it 
provides information on collisional, collective and electromagnetic processes in astro- 
physical and laboratory plasmas [13] and can help in the solution of various practical 
problems [14]. 

Since late seventies there have been several attempts to develop a theoretical model for 
understanding laser-matter interaction in ambient atmosphere but every model has 
some limitations. Knight [15] developed a theoretical model which was generalized 
by Bellantone et al [10,11] but its usage was restricted to relatively low intensity 
regime. Aden et al [12] developed a three dimensional model which matched well 
with their experimental data. Macfarlane et al [16] and Giuliani et al [17] had giver 
two separate models for the same experiment which was valid above 1 Torr but failec 
at higher pressure. Ananin et al [181 had given a model for plasma expansion at low 



works well at higher pressure. However, model is true for an instantaneous point 
explosion and with negligible radiation losses. Further, the model is insufficient for 
getting the detailed thermodynamic and hydrodynamic properties of the blast wave. 

The present paper deals with the case of intermediate intensity regime, 10 8 -10 12 
W/cm 2 and higher pressures. We take help of the modified version of the laser fusion 
code MEDUSA [20], MED101. This is an one dimensional Lagrangian code which 
calculates hydrodynamic and thermodynamic behaviour of the laser produced plasma in 
vacuum. MED101 is a well developed code which can take care of double and triple 
layered targets and can calculate X-ray laser gain. We have made modifications in this 
code to implement it for ambient gas. The Rankine-Hugoniot relations for shock waves 
are used to calculate the changing boundary conditions, that is temperature and pressure. 
The Navier Stokes equation is changed due to Reynolds stress arising due to turbulent 
boundary layer mixing and at a later stage due to diffusion and interparticle collisions. 
The thermodynamics of the system is also changed accordingly. The equations are solved 
using the finite difference method. 

The results simulated using modified MED 101 are compared with the experimental 
results of Puretzky et al [6] where the graphite target was irradiated by focussed beam of 
KrF laser and the plume expanded into a background of He, Ne, Ar and Xe gas. 



2. The numerical model 

2.1 In vacuum 

Here, we present a brief outline of the one dimensional Lagrangian code MEDUSA. A 
detailed description and the mathematics of the code is available elsewhere [20]. The 
plasma is assumed to consist of a charge neutral mixture of electrons and various species 
of ions. MEDUSA treats the thermodynamical properties of the electrons as one 
subsystem with internal energy u e per unit total mass, temperature T e , pressure p e and 
specific heat ratio % and so on, and the ions as second subsystem with a corresponding 
set of variables. Electric and magnetic fields are neglected and the two subsystems are 
coupled together by a common velocity and by energy exchange due to collisions 
between them. This is the one fluid two temperature model, the whole of the plasma is 
considered as one fluid in which electrons and ions have different temperature. According 
to the first law of thermodynamics 

d> dU dW 

^ __ I __ f-i \ 

dt dr dr ' (L) 
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The source terms Si and 5 e for ions and electrons respectively are written as 



(3) 
(4) 

Here H represents the flow of heat due to thermal conduction; K is the rate of energy 
exchange between ions and electrons; Y is the rate of thermonuclear energy release; J is 
the Bremstrahlung emission; X is the rate of absorption of laser light; Q s is the rate of 
viscous shock heating. 

Since plasma is considered to be a fluid, the motion of the plasma is governed by the 
Navier-Stokes equation 

P%=-VP> (5) 

where u is the velocity of plasma which defines the motion of the lagrangian coordinate, p 
is the total pressure due to electrons and ions together with artificial viscous pressure [20]. 
It is necessary to specify both hydrodynamic and thermodynamic boundary conditions 
at the inner and outer boundaries. For inner boundary u(r = 0) =0. For the outer 
boundary which is moving at a point r = RQ there can be various cases like 



2. (/?) = 0, 7i(K ) = 7i(0, T t (Ro) = r e (r), 

3.p(R Q }=p(t], 

4. (/?(,) = (*) (6) 

Boundary condition 1 is being the case of plasma expanding in vacuum. 

The equation of state assumes perfect gas equation of state for ions and Thomas-Fermi 
equation of state (with modified expression to give correct solid density) for electrons. 
These equations are solved along with the mass continuity equation by finite difference 
method to get temporal and spatial profile of plasma parameters viz. p, u, T\, T e , n\, n e and 
also the ionization state of the plasma. 

The original MEDUSA code is modified to MED 101 to calculate X-ray laser gain. The 
modified version includes atomic physics and radiation packages which can handle high 
and low Z materials, and a corrected equation of state for the electrons. The atomic 
physics package with the help of steady state CORONA model incorporates variable 
ionization in the plasma assuming equilibrium between collisional ionization and radiative 
recombination and depends only on electron temperature. It also includes time dependent 
form of the three body recombination which depends on electron density. In the older 
version of the code thermal electrons and ions were considered to be the only means of 
energy transport. The radiation package takes care of the enhanced energy transport due 
to suprathermal (fast) electrons and photons created at the critical density, whose mean 
free path is much longer than that of the thermal electrons resulting in fast energy 
transport [21]. The two equation of state available before were (1) perfect gas equation of 
state and (2) Fermi Dirac equation of state. Neither of them accounts for changes in 
ionization and hence the number of free particles remained constant. The Thomas-Fermi 
equation of state is introduced but it cannot deal with solid or liquid phase because it does 
not include binding forces. 'Corrected Thomas-Fermi' equation is introduced but it fails 
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Figure 1. Region I is expanding plasma, region II is compressed gas and region in is 
the undisturbed gas. 

governed by a perfect gas equation of state whereas for electrons one can choose any of 
the five states stated above. 

2.2 In ambient gas 

As the plasma moves at a velocity of 10 6 -10 7 cm/sec into the ambient gas it acts as a 
supersonic piston giving rise to a shock wave which moves at a faster speed than the 
piston [22]. If the gas pressure is high (say, > 20Torr) then there is a constant 
discontinuity between the piston and the swept up gas and there is hardly any intermixing 
of the two. A rough sketch of the process is given in figure 1 . Region I is the plasma 
region, region II is the compressed gas ahead of the plasma front and region III is the 
undisturbed gas where pressure.and temperature are PQ and TQ. Between region I and II is 
the plasma front which is referred to as supersonic piston. Between region II and III is the 
shock front which is moving faster than the plasma front. 

The velocity of the plasma front V p and the velocity of shock front V s are related by the 
expression [22] 



(7+1) 
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where a\ = 
the gas. 
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^RT is the velocity of sound in the medium and 7 the specific heat ratio of 



2.2.1 Boundary conditions. The boundary condition on the region I is defined by the 
pressure and temperature of the region n. Rankine-Hugoniot relations are used to 
calculate the temperature and pressure of region II. The situation is similar to figure 2(a). 
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Figure 2. (a) The moving shock front as seen from reference frame fixed to earth, 
(b) The stationary shock front as seen from reference frame moving with the shock 
front. 



Here the pressure (P c ) and temperature (T c ) in compressed gas region II across a 
stationary shock front are given by Rankine Hugoniot relations [22] 



where 7 is the specific heat ratio of the gas and M s = V s /a is the Mach number i.e the 
boundary condition is changing with time, the changed boundary conditions are 

P(R Q ) = P Q if M s < 1 
= P C ifM. s >l, 

T Q (Ro) = T Q if M s < 1 

= T C ifM v >l (10) 

If the Mach number is < 1 then shock waves are not formed and the boundary conditions 
are defined by PQ and TQ but if Mach number is > 1 then shock waves are formed and 
pressure and temperature are given by P c and T c respectively. 

Here the boundary condition for the temperature is applicable to the electron 
temperature only. It is because of the higher concentration of electrons at the plasma front 
due to their higher mobility than that of ions. Here the velocity is being calculated at 
every step. That is, after every small time step the code calculates the plasma front 
velocity and simultaneously it calculates the shock front velocity and then the 
corresponding pressure at the boundary and 'with that pressure it calculates the next 
plasma front velocity and this process continues till the end. 

2.2.2 Viscous drag. When plasma moves into an ambient gas it is almost similar to the 
case of hot jet coming out into ambient atmosphere. The moving plasma undergoes 
turbulent boundary layer mixing due to which it experiences viscous force on its 
boundary. 

According to Prandtl's mixing length theory the total Reynolds stress due to turbulent 
mixing in the boundary layer is given by the expression [24] 
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where e is kinematic viscosity, p is the density of the external gas, w max is the velocity of 
the center of flow and M m j n is at the boundary. x\ is an empirical constant to be determined 
from experiment, b is the width of the mixing zone. 

At the outermost boundary, velocity is zero. So M m in = and writing M TOax = u, 
eq. (11) simplifies to 

(12) 
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81 is the length of a particular mesh varying with time. The Navier-Stokes equation 
(eq. (5)) changes to 

P~ = -V/>-r(AO. (14) 

2.2.3 Diffusive force. As the plasma front expands into the ambient gas the pressure just 
inside the boundary decreases and the pressure outside increases. At a later time the 
pressure inside the plasma front is comparable to the pressure outside. At that time the 
piston becomes leaky and outside gas diffuses into the plasma which gives rise to 
collisions due to which the expanding plasma experiences a large resistive force 

C\ i 
PexA 
P 
PmlJ 

or 



Pext is the pressure just outside the plasma front and pi nl is the pressure just inside the 
plasma front, p is the density of the external gas. On addition of a diffusive term the 
modified Navier-Stokes equation reads 

dM T7 fAA P Y7 2 r(P\ n Ufi\ 

p = vp r(A/) Fdiff = Vp - px\u K\ p. (16) 



This equation contains three constants viz. *i, K and n to be determined from the 
experiments. 

3. Results 

With the modifications as above the simulation is carried out with an initial step of 
t = 10~ 18 sec having mesh size of 60. One of the simulated curves is compared with one 
of the experimental curves (say Xe) and constants x\ t K and n in eq. (16) are determined. 
These constants once determined are then used to get other simulated curves (in case of 
He, Ne and Ar) which are compared with respective experimental curves. The simulated 
curves are matched with the experimental data of Puretzky et al. Here graphite target is 
ablated by KrF laser (A = 248 nm, 30-ns FWHM, 2x2 mm 2 spot size on the target, 
fluence = 20 J/cm 2 ). The ablation is carried out in four different buffer gases at a 
pressure of SOOTorr. The temporal evolution of the plasma front as it propagates in 
He, Ne, Ar and Xe gas is shown in figures 3(a)-(d). The simulated results shown in 
solid line match well with the experimental results, shown in dashed line till 300nsec. 
The front propagates faster in lighter gas (like He) than denser gas (like Xe). It is 
interesting to note that after 300nsec the simulated value is more than that of the 




100 200 300 400 500 
Time(nsec) 

(a) 



100 200 300 400 500 
Time(nsec) 

(b) 




100 200 300 400 500 
Time(nsec) 

(c) 



100 200 300 400 500 
Time(nsec) 

(d) 



Figure 3. The position of .the plasma front versus time in four different ambient 
gases (a) He (b) Ne (c) Ar and (d) Xe. Experiment ( ), simulation ( ). 
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Figure 4. The temporal profile of (a) velocity (b) shock intensity of the plasma front 
in Xe gas as obtained from simulation. 



experimental one. This is because after 300 ns the diffusion increases to such an extent 
that a reflected shock is generated within the plasma plume which moves in the opposite 
direction, towards the target surface, damping the motion of the plasma as observed in 
the experiment. If there is mixture of fluids whose composition is different at different 
points then the equations of fluid dynamics have to be considerably modified [25]. So a 
different set of equations are needed to model this part which is beyond the scope of the 
paper. 

In figure 4(a) the temporal profile of plasma front velocity with Xe as buffer gas as 
predicted by this numerical model is shown. The velocity decreases with time but at 
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Figure 5. Figure 3(d) and figure 4(a) is replotted with a common time axis. 



shock wave going into the plasma as a function of time as observed in simulation. The 
shock intensity is maximum at around 290ns, that is at the same time when velocity 
suddenly drops and increases again. Experimentally it was observed that a shock wave 
had gone into the plasma front at around 300 ns causing instability in the velocity of the 
plasma front. The simulated result match very well with the experiment. It has also been 
observed in simulation that the intensity of shock wave going into the plasma in case of 
He is orders of magnitude lower than that of Xe. It is because Xe being a denser gas 
gives a greater thrust to the moving plasma compared to that of He. Figure 3(d), that 
is position of the plasma front expanding in Xe and figure 4(a) that is the velocity of 
the plasma front in the same ambient condition is replotted in figure 5, with a common 
time axis, position and velocity axis are shown along the Y axis in left and right 
respectively. The whole graph is divided into seven different regions to show that slope of 
the position graph obtained by simulation gives the simulated velocity graph. For e.g. in 
region (f) the slope of the position curve suddenly tends to 0, the velocity also tends to 
at that point as seen in the velocity curve where as in region (d) and (g), the velocity is 
constant. 

The simulated temporal profile of temperature, pressure, density and electron density 
of the plasma front in case of Xe (dashed line) and He (solid line) as buffer gas is shown 
in figures 6(a)-(d). Xe being denser than He gives higher resistance to the moving plasma 
resulting in higher temperature, pressure, density and electron density of the plasma front. 
The temperature of the plasma front is highest at around 30 nsec (figure 6(a)) which is the 
pulse duration time, after which it decreases due to the cooling effect of the inert gases. 
The density profile in case of Xe (figure 6(c)) show a sudden rise at around 290 nsec due 
to shock wave coming into the plasma . Theoretically, according to perfect gas equation 
of state, rise in density at around 290 nsec (figure 6(c)) and fall in its temperature at the 
same time (figure 6(a)) gives no significant change in pressure. This is what is observed 
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Figure 7. The simulated position of the plasma front as obtained from (a) eq. (5) i.e. 
no change in Navier-Stokes equation, (b) eq. (14) i.e. on adding viscous drag term 
only, (c) eq. (16) i.e. on adding viscous drag term and diffusive force term. 



in the simulated pressure profile at 290ns in figure 6(b). The inset in figure 6(c) and 
figure 6(d) show the density profiles and electron density profiles for He in a magnified 
scale. It shows that density as well as electron density increases similar to that of Xe but 
increase is gradual and low in case of He compared to Xe where change is abrupt and 
high because of its higher density. 

Figure 7 shows the effect of different terms which are added in the Navier Stokes 
equation other than the boundary conditions. The viscous drag term plays the major role 
in decelerating the plasma front. The diffusive term comes into play at around 300 ns and 
is much less compared to viscous drag. 

The model successfully predicted the results of another experiment by Gu et al [5] also. 



4. Conclusion 

The laser matter interaction in the presence of an ambient gas has been modeled 
numerically by making appropriate modifications in the one dimensional laser fusion 
Lagrangian code MEDUSA. The model is applicable for the high pressure of ambient gas 
where diffusion and interparticle collisions are not dominant in the beginning of the 
plasma expansion. The simulated result such as position of the plasma front, instability in 
velocity and the time of evolution of shock wave going into the plasma front are in 
reasonable agreement with the experimental observations [6]. Using the modified code the 
model developed can be used to calculate spatial and temporal profile of the plasma 



References 

S Lybchenko, I K Navikov and S P 



Fnln ' 

Frolov, Sov. J. Quantum Electronics 21, 787 (1991) 

[2] M A Liberman and A L Velikovich, Phys. Fluids Bl, 1271 (1989) 
[3] H S Kim and A S Kwok, Appl. Phys. Lett. 61, 2234 (1992) 
[4] David Geohegan, Appl. Phys. Lett. 60, 2732 (1992) 

C Chen ' Z Y Wan S and P K Lin, Appl Phys. B58, 143 



[6] A A Puretzsky, D B Geohegan, R E Hauflerm, R L Hettich, X Y Zheng and R N Compton in 

AIP Conference Proceedings 288, Laser Ablation: Mechanisms and Applications-lifted by 

J C Miller and D B Geohegan (AIP Press, New York, 1994) pp. 365 
[7] E Hontzopoulos, D Charalambidis, C Fortakis, Gy Farkas, Gy Horvath and Cs Toth Opt 

Commun. 67, 124 (1988) 

[8] Daniel Fried and Gene P Reck, J. Phys. 24, 1065 (1991) 
[9] T Pisarczyk, A Farynski, H Fiedorowicz, P Gogolewski, M Kusnierz, J Makowski, 

R Miklaszewski, M Mroczkowski, P Parys and M Szczurek, Laser and Part. Beams 10 

767 (1992) 

[10] R Bellantone and Y Hahn, J. Appl. Phys. 76, 1436 (1994) 
[11] R Bellantone and Y Hahn, J. Appl. Phys. 76, 1447 (1994) 
[12] M Aden, E Beyer, G Herziger and H Kunze, J. Phys. D2S, 57 (1992) 
[13] B H Ripin, C K Manka, T A Peyser, E A Mclean, J A Stamper, A N Mostovych, J A Stamper, 

J Grun, K Kearney, J R Crawford and J D Huba, Laser and Part. Beams 8, 183 (1990) 
[14] D L Pappas, K L Saenger, J Bruley, W Krakow, T Gu, R W Collins and J J Cuomo, J. Appl. 

Phys. 71, 5675 (1992) 
[15] C J Knight, AIAA J. 17, 519 (1979) 

[16] J J Macfarlane, G A Moses and R R Peterson, Phys. Fluids B17, 635 (1989) 
[17] John L Giuliani, Jr. and Margaret Mulbrandon, Phys. Fluids B17, 1463 (1989) 
[18] OB Anan'in, Yu A Bykovskii, Yu V Eremin, A A Zhuravlev, O S Lyubchenko, I K Novikov 

and S P Frolov, Laser Phys. 2, 711 (1992) 
[19] L I Sedov, Similarity and dimensional methods in mechanics (Academic, New York, 1959) 

G I Taylor, Proc. R. Soc. London. A201, 175 (1950) 

[20] J P Christiansen, DETF Ashby and K V Roberts, Comp. Phys. Commun. 7, 271 (1974) 
[21] T D Beynon and E H Smith, Proc. 2nd International Conference on radiative properties of hot 

dense plasma (Sarasota, Florida, 1983) 

[22] James E A John, Gas dynamics, second edition, (Allwyn and Bacon, INC, Boston, 1984) 
[23] Zhaoyuan Han and Xiezhen, Shock dynamics (Science Press, Beijing Kluwer Academic 

Publishers, 1993) 

[24] H Slitching, Boundary layer theory (Pergamon Press, London, 1955) 
[25] Landau and Lifshitz, Fluid mechanics (Pergamon Press, New York,1984) 



Pramana - J. Phys., Vol. 50, No. 1, January 1998 73 



Absence of neutron emission during interaction of 
deuterium with metal at low energies 

A SHYAM and T C KAUSHIK 

Neutron Physics Division, Bhabha Atomic Research Centre, Mumbai 400 085, India 
Email: npd@magnum.barctl .ernet.in 

MS received 13 June 1997; revised 18 October 1997 

Abstract. Technique and instrumentation to detect reliably, multiplicity of neutrons emitted in 
sharp bursts (< 100 us) has been developed where a burst of as low as 15 neutrons and continuous 
emission of ~ 10" 1 neutron/s may be detected. Using this technique, attempts were made to detect 
neutron emission from various experiments in which anomalous nuclear effects (or what is 
commonly referred to as cold fusion) may be expected to occur. No neutrons, above our detection 
threshold, were detected in the recent series of experiments. 

Keywords. Anomalous nuclear effects; cold fusion; deuterated solids; burst neutrons; pulsed 
neutrons; neutron detection. 
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1. Introduction 

Several investigators, including the present authors [1,2], have reported on the emission 
of neutrons when deuterium is loaded into palladium or titanium by electrolysis or by 
other techniques [3-16]. These experiments are usually referred to as anomalous nuclear 
effects in deuterium-solid systems or cold fusion. The neutrons were reported to be 
emitted continuously, as well as in sharp bursts presumably by the D-D fusion reaction 

" 2 D] + 2 Di -> 3 He 2 + ni (2.45 MeV), (1) 

where 2 D] , 3 He 2 and % are deuteron, helium-3 and neutron respectively. 

2. Neutron detection systems 

In order to confirm and understand the nature of observed neutron emissions, simple, 
relatively inexpensive and high sensitivity systems had been developed [17] in two 
configurations, to characterize the burst and continuous neutron emissions. These 
detection systems were designed and optimized based on neutronics numerical 
simulations [18] and then through experiments. Hydrogenous materials (epoxy and 
polyethylene) were used to thermalize, back scatter and stretch (~ 100 us) the neutron 
burst. The thickness of hydrogenous material (to thermalize the neutrons), between 
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Figure 1. Schematic of the first neutron detecting system. 



material was at the back-side of the detectors. This material helped in back-scattering the 
escaped (signal) neutrons into the detectors. It also, partially, shielded the detectors from 
background neutrons. 

In the first system developed here, 12 boron-tri-fluoride (BFs) proportional detectors 
were arranged into three independent channels (see figure 1). Two detector arrays were 
used to record signal and the third to monitor background. The neutrons were 
continuously detected in 100|j.s and 5s interval in order to detect burst as well as 
continuous emissions. The second system had 16 BFa detectors in two concentric rings 
(see figure 1). In this configuration, the second ring acted both as background and signal 
detector. In both the systems weather-proof connectors were used. The preamplifiers were 
electromagnetically shielded and kept close to detectors. The system was reliable enough 
to work satisfactorily even in 100% condensing humidity (of Mumbai monsoon). 

Each channel had separate counting electronics (though powered from single source 
and grounded at same point). The counts from each channel were recorded in a personal 
computer controlled scalar, communicating to an IBM PC/486 compatible computer 
through RS-232 serial interface. 

The efficiencies were measured using a standard californium-252 source (average 
neutron energy ~2.2MeV, close to D-D fusion neutron energy). The efficiencies of the 
channels, in the first system were ~3.9 0.4%, 4.6 0.4% and 3.9 0.4%, when source 
was kept close to the detector channels. When the source was kept in the experimental 
space the efficiencies were ~3.5 0.4%, 4.2 0.4% and 0.5 0.4% for the two signals 
and the background detector banks, respectively. The average background/s in the three 




Figure 2. Schematic of the second neutron detecting system. 

detector systems were 0.062 0.005, 0.076 0.006 and 0.063 0.005 counts per 
second (cps). For the second system, the efficiencies and the backgrounds of the three 
channels (with calibration source kept in the central experimental area) were 5.1 0.4%, 
4.9 0.4%, 3.0 0.4% (outer ring), 0.060 0.005, 0.063 0.005, 0.067 0.006 cps 
respectively. However, when the calibration source was kept at periphery of the detecting 
system, at six different locations, all around the circumference, the average efficiencies 
were 1.1 0.3%, 1.0 0.3% and 2.5 0.4% (outer ring). Therefore in this setup, the 
outer ring of detectors, served the dual purpose of signal as well as background channel. 



3. Statistics of detection 

3.1 Continuous emission 

If the counting is done for a very long duration (as in cold fusion experiments), in other 
words if there is a large number of counting intervals, there is always a probability of 
getting counts in an interval which are several a above the background. Therefore the 
data was analysed by comparing the experimentally observed frequency distribution with 
Poisson distribution 

P p (c}=e- a -a c lc\, (2) 

where P p is the frequency, predicted by the Poisson distribution, of getting counts c in a 
counting time-interval, a being mean of counts. The goodness of fit between observed 
and Poisson frequency distribution was determined by x 2 (chi-square) test. 

As we have used more than one detecting channels, the binomial probability (for 
isotropic emission) of getting M counts in detecting channel 1 and N counts in the 
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This probability is normalized to give Pt, = 1 for both the counters giving equal counts 
(M = N}. This indicates that probability of getting more than twice (or half) the counts in 
one channel (N) as compared to the counts in another channel (M) is insignificant, unless 
the counts recorded are very small (M < 2). By this method the spurious counts, if 
generated in one channel, due to some defect, could be ignored. The thumb rule criteria 
used to look for genuine neutron emission was 

(i) the counts in one of the signal channel should not be greater than twice (or less than 
half) the other channel (after normalizing for channel efficiencies), 

(ii) the counts in the background channels (a) should be close to background counts, in 
the first system and (b) should be much less than counts in the signal channels, in the 
second system, 

(iii) the frequency of the event(s) should be higher than predicted by Poisson distribution. 

3.2 Burst emission 

For burst neutron emission, if the average background per second is b, and r is the 
counting interval, then the Poisson probability P s of getting j counts in the interval is 

To get multiplicity of at least 1 (j = 2), from continuous emission (in the counting 
interval of 100 jis and background of ~0.07/s - as encountered in our experiments) the 
probability is 2.2.10" 11 or 1 multiplicity count should be observed every 2 months. 
Therefore observing even a single multiplicity count indicates a burst emission. However 
as the cosmic ray background also generates bursts, the criteria to identify neutron bursts 
from the experiments, should be similar to conditions stated for continuous emissions 
(described later). 



4. Experiments and results 

Experimental configurations similar to those which were previously performed, and had 
reported neutron emissions [5,8,10,13,19-24] were repeated further. The neutron 
emissions reported in these experiments, were high enough to be detected easily, by the 
present detecting systems. 

One of the most extensively investigated devices is being described in detail here. The 
device studied was commercially available palladium(Pd)-nickel(Ni) electrolysis cell (or 
the Milton Roy cell as it was procured from M/s Milton Roy, Ireland). Other reasons 
for choosing this configuration of electrolytic cell as an illustrative example are (a) it was 
reported [25] to have produced large neutron emissions in continuous as well as in burst 



third run, there appeared to be a low, but definite increase in burst and continuous counts 
during electrolysis but in a subsequent ran counts were observed to be lower during the 
electrolysis. 

4. 1 Commercial palladium-nickel electrolyser 

The commercially available Milton Roy cell had 16 Pd tubes (3 x 10~ 3 m outer 
diameter, 0.2 m long), coated with Pd black and arranged in a circular array, as cathode. 
The total surface area of the tubes was ~ 3 x 10~ 2 m 2 . The Pd cathode was surrounded 
on both sides by two concentric, cylindrical Ni anodes. Three cells of this type were 
operated for a total five times. These experimental runs (or series) have been labelled as 
MR1, MR2, MR3, MR4 and MRS. 

4.2 MR1 series of experiments 

The cell was operated with electrolyte consisting of 5 molar (M) sodium deuteroxide 
(NaOD) dissolved in heavy water (D2O). The current was maintained at 60 A 
(2000 A/m 2 ). The results obtained for MR1 series have already been reported [9, 10]. 
However it may be noted that during this series all the counters were configured into one 
segment, connectors were not weather-proof and the preamplifiers did not have 
electromagnetic shielding. 

4.3 MR2 series of experiments 

This series of experiments were replication of the MR1 series. The cell was again 
operated at 60 A current, using 5 M NaOD solution in D2<3 as electrolyte . However in 
this series, the BF 3 detectors were segmented to form three independent channels as 
shown in figure 1. Along with this system, proton recoil liquid scintillator detectors, with 
pulse shape discrimination capability were also used. The latter could detect fast neutrons 
(energy threshold ~!MeV) with a detection efficiency of ~40% till 14MeV neutron 
energy. With these detectors it was possible to detect fast, continuously emitted, neutrons 
with a detection threshold of 0.7 n/minute [26]. 

During this series of experiments, no neutron emission could be observed, either in 
continuous [26] or in burst mode. It may be noted that this (Milton Roy) cell was 
equivalent to at least 16 simple cells [3]. Therefore the neutron yield expected from MR 
cell was at least one order of magnitude higher than reported from a simple cell [3] and 
thus well above the detection threshold of our detecting system. 

4.4 MRS series of experiments 

During the third series of experiments (and subsequent runs), the second type of BFa 
neutron detecting system (see figure 2) was employed. As other investigators had 
reported to observe neutron only after several weeks of loading and using lithium 
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system's (without cell) background was observed. In the second phase of experiments, the 
cell was operated at 10 A for the next 30 days, using saturated solution of LiOD in D2O as 
electrolyte. In the third phase, the cell was run for (next) 15 days using saturated solution 
lithium hydroxide (LiOH) in natural water (H 2 0) as electrolyte. 

The burst and continuous counts per day are plotted in figures 3a, b respectively. It 
appears from these plots that neutrons were getting produced (burst as well as 
continuous), not only during D 2 O electrolysis (as compared to blank run) but also, to a 
lesser extent, during H2O electrolysis [2]. 

4.5 MR4 series of experiments 

The fourth set of experiments (MR4 series) were performed in an attempt to reproduce 
the results of MR3 series with more controls. In this series, for the first ten days data was 
recorded without cell, for the next ten days data was obtained with cell filled with 
saturated solution of LiOH in H 2 O, but without electrolysis. For next 10 days data was 
recorded with cell filled with saturated solution of LiOD in D 2 O, but without electrolysis. 
For the last ten days, the cell was operated at 10 A using saturated LiOD in D 2 O as 
electrolyte. As can be seen from figure 3d, the continuous counts observed were lowest, 
during electrolysis (carried out at 10 A) using LiOD-D2O electrolyte. This was contrary 
to the behaviour observed during MR3 series. 

4.6 MRS series of experiments 

Due to the contradictory results obtained in MR3 and MR4, long duration experiments, 
MRS series of experiments were done for shorter durations. Following experiments were 
performed for 24 h each: (a) blank, (b) cell only, (c) cell filled with H 2 O, (d) cell filled with 
saturated solution of LiOH in H 2 O, (e) electrolysis of this cell at 10 A, (f) cell filled 
with D 2 O, (g) cell filled with saturated solution of LiOD in D 2 O, (h) electrolysis of this 
cell at 2 A, (i) electrolysis of the same cell at 4 A, (j) electrolysis of the same cell at 8 A, 
(k) electrolysis of the same cell at 16 A and (1) electrolysis of the same cell at 32 A. 

In this set of runs, no neutrons, in bursts or in continuous mode, could be observed 
above the detection threshold of the counting system. All the counts were within the 
statistical variations of the background (the statistical variations in this series of 
experiments were higher because of shorter duration of the data collection). 

5. Discussion 

The results obtained in MR4 series were contradictory to those obtained in MR3 series. 
Such observations may however be explained if the neutrons are assumed to have been 
produced by cosmic rays (secondary) impinging on the cell and producing neutrons by 
(nuclear) reactions with Pd, Ni and lithium in H 2 O-LiOH cell. Moreover in D 2 O-LiOD 
cells, additional number of neutrons may be produced by disintegration of deuterons by 
cosmic ravs. As the cosmic rav intensity varies oeribdicallv 121], it may be conjectured 



counts were lowest during D 2 O-LiOD electrolysis. Since the neutrons were produced by 
cosmic rays impinging on the cell, the signal channels recorded more neutrons than the 
background channel. 

The number of neutron burst counts were higher, in the case of D 2 O-LiOD cell both in 
MR3 and MR4 runs. However in the MR4 run, such counts were high even before 
electrolysis and remained approximately the same afterwards (see figure 3c). Hence these 
bursts could not be attributed to electrolytically induced nuclear reactions. Such bursts 
were perhaps again produced by large (high multiplicity) cosmic ray showers. More 
number of bursts observed with the DiO-LiOD cell may be accounted through possible 
extra neutrons produced by deuteron disintegration. Muons are also known to be 
generated in cosmic showers which may further catalyse fusion reactions and contribute 
[28] to the registered low level neutron yield in the LiOD-D2O cells. 

6. Conclusion 

Reliable and high sensitivity detecting systems were used to search for burst as well as 
continuous neutron emission from electrolysis of D2O-LiOD. No conclusive evidence of 
neutron emission (pulsed or continuous) could be obtained from a number of such runs. 
A neutron detection system needs to be properly designed to avoid malfunctions which 
may give rise to misleading signal of neutron emissions from an experimental source. In 
the long duration experiments, the change in cosmic ray background, also appears to be 
an important effect which could reasonably explain many contradictory observations 
from such experiments. 
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Abstract. A toroidal grating monochromator (TGM) based beamline was designed earlier and 
it is being^ setup on the 450 MeV Indus-1 synchrotron radiation source. This beamline operates in 
40- 1000 A wavelength range and gives a moderate resolution (~500) at high throughput. Now-a- 
days, spherical gratings are giving very high resolution in some complex monochromators. We have 
evaluated the use of spherical gratings in the simple rotation mount of a TGM. The motivation has 
been to evaluate if the existing TGM beamline, without any hardware changes, can be used with the 
much cheaper spherical gratings. The two monochromators, employing toroidal and spherical 
gratings respectively, have been compared by performing detailed ray tracing calculations. Our 
study shows that, with the same mechanical hardware, a real spherical grating monochromator gives 
better performance compared to the conventional TGM with respect to the spectral resolution and 
the flux. Furthermore, for moderate resolution this improvement is realized even when the exit slit 
is not moved to keep the spherical grating in-focus. 

Keywords. Synchrotron radiation; beamline; monochromator; soft x-ray. 
PACS No. 42.10 



1. Introduction 

A number of synchrotron radiation facilities exist world wide, many of which emit 
radiation in the soft x-ray and vacuum ultra violet (XUV) wavelength range (10-1200 A). 
In this wavelength range, grazing incidence optics in ultra high vacuum (UHV) 
environment has to be used. Many designs of the monochromators to be put up on beam 
lines for this wavelength range exist [1,2] and are based on plane, toroidal or spherical 
gratings. Many of these designs require stringent closed coupled movements of the 
optical elements in the UHV environment and are therefore difficult to realize in practice 
and are quite expensive. In most cases, these complicated designs are essential to derive 
the best resolution from the monochromators. However, there are many application areas 
like photoelectron spectroscopy, reflectometry and radiometry where the requirement is 
not for very high resolution but of good throughput at moderate resolution (A/AA ~500- 
1000). Relatively simple monochromator designs can deliver these requirements leading 
to a great deal of simplicity in the instrumentation. One such design [3] that has been 
quite popular during the eighties is the use of toroidal gratings in the constant deviation 
geometry (TGM) wherein a simple sine bar drive rotates the grating about its axis for 
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figure errors of ~5 arcsec. These tangent slope errors severely limit the resolution 
achievable with these monochromators. However, since there is only one optical element 
between the entrance and exit slits, these monochromators yield high throughput. 

On the other hand, spherical and plane surfaces are easier to fabricate with better 
figures. Spherical surfaces with 1 arcsec slope error are quite common. However, earlier 
not many monochromators were designed using spherical gratings in the XUV 
wavelength range because spherical surfaces suffer from strong astigmatism and coma 
aberrations at grazing incidence. In the last decade, a few good designs [4-7] for 
spherical grating based monochromators have come up which give very high spectral 
resolution. These designs generally require complicated mechanical movements and very 
small slit widths. In some designs [6, 7], the exit slit, and also at times the entrance slit, is 
required to be moved to satisfy the Rowland circle condition. This requirement of moving 
the slits with wavelength scan is, at times, not desirable for certain experiments. 

A beamline based on a TGM has been designed [8] and is being set up on the Indian 
synchrotron facility Indus- 1 [9] that is coming up in our centre. This beamline operates in 
40-1000 A wavelength range and gives moderate resolution (~500) at high throughput. 
Since now-a-days, spherical gratings are giving very high resolution in Dragon [5] and 
focusing-spherical grating type of monochromators [7], we have evaluated the perfor- 
mance of spherical gratings in the simple rotation mount of a TGM. The motivation has 
been to evaluate if the existing TGM beamline, without any hardware changes, can be 
used with the much cheaper spherical gratings in a later upgrade of the beamline. To 
evaluate the feasibility of such a scheme, we have evaluated the configuration for the 
monochromator in which toroidal gratings in the TGM mount are replaced by spherical 
gratings, and suitable preoptics [10] is used which decouples the meridional and sagittal 
focus positions to eliminate the astigmatism of the spherical grating. Using ray tracing 
calculations, the performance of this spherical grating based monochromator is evaluated 
and compared with the conventional TGM monochromator designed earlier. The details 
of this intercomparison is given in this paper which shows that the monochromator with 
spherical gratings give better performance both in terms of resolution and throughput. 
Moreover, though this intercomparison has been done for the Indus- 1 synchrotron source, 
the conclusions drawn are fairly general and will hold for any other XUV synchrotron 
source also. 



2. Optical configuration 

The requirement put on the beamline is that it (i) produces monochromatic (A/AA ~500- 
2000) radiation in 40-1000 A wavelength range, (2) provide a relatively high intensity 
beam of small size at the sample position, and (3) yield a spatially fixed beam after the 
monochromator. 

The detailed optical design of the beamline using TGM was derived earlier [8]. 
Parameters for the SGM, to fit in the same mechanical hardware, have been determined 
by a semi-empirical optimization procedure. It involves the optimization algorithm of 
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resolution (AA), integrated over the desired wavelength range, was minimized. However, 
on comparing this resolution with that of TGM it turned out that, though the overall 
integrated resolution over the whole wavelength range was better than TGM, the resolu- 
tion at higher energies (i.e., lower photon wavelengths) was worse than the TGM. Hence 
the SGM was reoptimized semi-emperically so as to get similar or better resolution at 
each wavelength in the wavelength range of interest. During the optimization, we were 
constrained to keep the entrance arm length, the exit arm length and the deviation angle 
fixed, in order to fit the SGM into the hardware of TGM, and only the radii of the grating 
was varied. A total of three gratings are used to cover the full wavelength range. Various 
parameters of the monochromator are given in table 1 . 

The optical layout of the two configurations is given in figure 1. Configuration 1 
comprises of a premirror MI, an entrance slit Si, a toroidal grating TG, an exit slit 82 and 
a post mirror M2- The prefocussing toroidal mirror is kept at a distance of 4500 mm from 
the source, which deflects the SR vertically by 9. The mirror demagnifies the SR source 
by a factor of 2 and focuses both in sagittal and meridional direction at the entrance slit. 
The grating in the TGM is vertically dispersing at a constant deviation of 162. A second 
toroidal mirror is used to refocus the exit slit image onto the sample. 

Configuration 2 also comprises of a premirror MI, an entrance slit Si, a spherical 
grating SG, an exit slit S 2 and a post mirror M2- The premirror MI reflects the beam 
vertically by 9, but unlike the premirror of configuration 1, it focuses the beam 
horizontally and vertically on different positions. The two radii of this toroidal mirror are 
so chosen that it focuses the beam horizontally on the exit slit and vertically on the 
entrance slit of the optimized monochromator. The ideal monochromator for this type of 
preoptics is a cylindrical grating [5] which does not modify the sagittal focus position of 
the premirror. Use of a spherical grating having a non-zero curvature in the sagittal 
direction brings in the horizontal focus point from the desired exit slit position. Since the 



Premirror Slit 

Ml SI 



Slit Post mirror 
S2 M2 



ELEVATION 




4500 


3210 
t 


3210 
< > 


1836 
> 



(a) 




4500 


3210 


3210 
( 


1836 
( > 



(b) 



All dimensions are in mm 



Figure 1. Optical layout of the two configurations (a) based on toroidal grating; (b) 
based on spherical grating. 

grating rotates about its axis for wavelength scan, this influence of the grating on the 
horizontal focus position is wavelength dependent. We have increased the horizontal 
radius of the toroidal premirror empirically so as to get the sagittal focus position at the 
exit slit at a wavelength in the middle of the wavelength range of interest. A toroidal post 
mirror is again used to refocus the exit slit image onto the sample. 

3. Ray tracing studies 

The two configurations for the beamline were ray traced using the ray tracing code RAY 
[11]. The radiation source is taken to be the dipole of the 450 MeV Indus-1 storage ring 
and its various parameters are given in table 2. The vertical divergence of the SR beam 
from a bending magnet source varies as the wavelength is scanned and is in the range of 
2-9 mrad for Indus-1. In the horizontal direction, the beamline accepts a fan of lOmrad. 
Various parameters of the gratings and mirrors given in table 1 are used for ray tracing. 
The tangent slope errors of 5 arcsec on toroidal grating and 1 arcsec on spherical grating 
have been considered which correspond to the typically achievable manufacturing 
tolerances for these grating shapes. 
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Figure 2. Spectral resolution (/ AE) as a function of energy for slit settings of 40, 
200 and 500 urn. 

The resolving power of the two monochromator configurations for slit settings of 40, 
200 and 500 urn, respectively is shown in figure 2. The figure shows that in the low- 
energy range of each grating, the two configurations give similar resolution while for the 
rest of the spectral ranee the SGM configuration sives much hieher snectral resolution 
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Figure 3. Spot diagrams at 149.4, 150.0 and 150.6 eV for (a) TGM configuration, 
and (b) SGM configuration. Also shown is the intensity distribution in the horizontal 
and vertical directions. 



than the TGM configuration. For instance, with slit size of 40 um, the SGM configuration 
gives a peak resolving power of 2300 as against 900 available from the TGM configura- 
tion. The peak in the resolution curve, with narrow slit widths, is sharper in the case of 
spherical gratings. This is a manifestation of the fact that the spherical gratings give very 
high resolution when the Rowland circle condition is satisfied (at the peak). This feature 
is exploited many times by the beamline designers by choosing the peaks, by appropriate 
choice of line densities and grating curvature, to lie in the near edge region of specific 
elements. Figures 3a and 3b show the spot diagrams at the exit slit planes of TGM and 
SGM respectively for photon energies 149.4, 150.0, 150.6 eV. While the TGM configura- 
tion is not able to resolve the three energies, with the SGM configuration these energies 
are seen to be well separated. 

RAY program also calculates the flux throughput of the beamline by considering the 
wavelength dependent radiation emission from the Indus- 1 bending magnet, the absolute 
reflectivities of the optical elements, figure errors on all optical elements including the 
gratings, and the geometrical loss in the beamline. The grating efficiency could not be 
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Figure 4. Photon flux at the sample position for the two configurations of the 
beamline. Slit settings have been kept at 200 jam. 

calculated, as the parameters of the commercial gratings employed were not available. In 
the soft x-ray regime, the grating efficiency is typically in 10 to 30% range and therefore 
for flux calculations we have assumed it to be 20% over the whole wavelength range. The 
flux available at the sample position from the two configurations is shown in figure 4 
from which it is clear that the flux throughput of the SGM configuration is higher than 
that from the TGM configuration. These ray tracing results therefore show that if the two 
configurations of the beam line are operated in exactly identical conditions then the SGM 
configuration gives better performance both in terms of resolution and flux. 



4. Discussion 

The principal aberration of a spherical grating in grazing incidence is the astigmatism, 
which gets compensated in configuration 2 by the use of a proper toroidal premirror 
which focuses sagitally and meridionally at separate points. The other major components 
of the aberration in grazing incidence geometry are the defocus, astigmatic coma and 
coma [13]. Through the optimization process, the monochromator is brought to focus at 
two wavelengths. Then at the in-between wavelengths, the monochromator remains 
approximately focused. The. defocus component of the aberration can be eliminated by 
moving the exit slit, as the wavelength is scanned, to keep the monochromator in-focus 
[13]. We have also determined analytically the various components contributing to the 
wavelength resolution and figure 5 shows the full width at zero height (fwzh) resulting 
from different terms for a SGM. This figure shows that, over most of the wavelength 
range, the astigmatic coma dominates the aberration limited resolution and, therefore, 
moving the exit slit to keep the monochromator in-focus will not improve the resolution 
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Figure 5. Various components of the aberration limited resolution (full width at zero 
height) derived analytically for the SGM configuration. 

significantly. Thus movement of the exit slit can be avoided which leads to a great deal of 
simplicity in the experimental station. Keeping the monochromator in-focus at all 
wavelength will help if the relative proportion of astigmatic coma is reduced compared to 
the defocus contribution. This can be done by masking the non-dispersive aperture of the 
grating since astigmatic coma depends on the square of the non-dispersive aperture. 
However, this will drastically reduce the throughput of the monochromator. 

5. Summary 

It is shown that the use of spherical gratings in a constant deviation TGM mount gives 
better performance if proper preoptics is used to reduce the astigmatism. Detailed com- 
parison of this monochromator with a conventional TGM is performed using ray tracing 
and is presented in this paper. If one is interested in a moderate resolution of ~1000 then 
one can get the improved performance even without moving the exit slit position. Use of 
spherical gratings have the obvious advantages - they are easier to fabricate to high 
accuracy and because of symmetrical shape they are easier to align in the beamline. 

Acknowledgement 

We thank BESSY, Berlin for kindly providing the source code of RAY program. 

References 

[1] R L Johnson, Nucl. Instrum. Methods A246, 303 (1986) 

[2] R L Johnson, in Handbook of synchrotron radiation edited by E E Koch (North-Holland, 

Amsterdam, 1983) vol. 1, p. 173 
[31 B P Tonner and E W Plummer, Nucl. Instrum. Methods 177, 153 (1980) 



|4J b U Brown, K Z Bachrach and Neil Lien, Nucl. Instrum. Methods 152, 73 (iy 78) 

[5] C T Chen, Mid. /nsrrwm. Methods, A256, 595 (1987) 

[6] H A Padmore, SPIE Proceedings 733, 253 (1986) 

[7] W B Peatman, J Bahrdt, F Eggenstein, G Reichardt and F Senf, Rev. Sci. Instrum. 66, 28 
(1995) 

[8] K J S Sawhney and R V Nandedkar, Pramana - J. Phys. 42, 49 (1994) 

[9] S S Ramamurthi and G Singh, Nucl. Instrum. Methods A3S9, 15 (1995) 
[10] W A Rense and T Violett, J. Opt. Soc. Am. 49, 139 (1959) 
[11] W R Mckinney and M R Howells, Nucl. Instrum. Methods 172, 149 (1980) 
[12] F Shaefers, RAY program, BESSY, Berlin (unpublished) 

[13] J B West and H A Pdmore, in Handbook of synchrotron radiation edited by G V Marr (Nor 
Holland, Amsterdam, 1987), vol. 2, p. 21 



journal of physics Vol. 50 (No. 1), January 1998 

Contents 



Universality in the effective pair interaction of d-shell metals - Compressibility 

and vacancy formation energy of 3d-liquid metals 

A K Karmakar and R N Joarder 1 

Evidence of polaronic charge carriers in oxide superconductors 

CM Srivastava 1 1 

Electrical properties of a-antimony selenide 

Sanjeev Gautam, D K Shukla, Shelly Jain and N Goyal 25 

Monte-Carlo optimization of correlated helium wave function 

Sabyasachi Kar and Puspajit Mandal 35 

Total electron scattering cross sections for nitrous oxide at low electron energies 
P Rawat, K P Subramanian and Vijay Kumar 45 

Hyperspherical calculation for electron hydrogen atom ionization collisions in 

a symmetric form J N Das 53 

Numerical modeling of laser matter interaction in the presence of an ambient gas 
A Neogi and R K Thareja 63 

Absence of neutron emission during interaction of deuterium with metal at low 
energies A Shyam and T C Kaushik 75 

Use of spherical gratings in a toroidal grating mount for vuv-soft x-ray syn- 
chrotron radiation K J S Sawhney, A I Munshi and R V Nandedkar 85 



Indexed in CURRENT CONTENTS ISSN 0304-4289 



